Method and system for part design using heterogeneous constraints

ABSTRACT

A method of classifying design criteria includes receiving design criteria for a product part. The criteria comprise one or both of performance and manufacturing criteria. The design criteria are sorted into different classes of one or both of one or more objective functions and one or more constraints based on when they can be satisfied or optimized. Constraint violations are determined. A design workflow is produced to generate one or more designs of a part to comply with one or more of satisfying constraints and optimizing objective functions.

RELATED APPLICATIONS

This application is a continuation of U.S. patent application Ser. No.16/561,633, filed Sep. 5, 2019, to which priority is claimed pursuant to35 U.S.C. § 119(e), and which is incorporated herein by reference in itsentirety.

STATEMENT REGARDING FEDERALLY SPONSORED RESEARCH AND DEVELOPMENT

This invention was made with Government support under DARPA contractHR0011-17-2-0030. The government has certain rights to this invention.

TECHNICAL FIELD

The present disclosure is directed to the design of mechanical parts.

BACKGROUND

Generating practical designs involves simultaneously reasoning aboutmanufacturing, assembly, materials, and physics. Computational services,herein called solvers, are used to assist this process by analyzing acandidate design to see if it meets some performance criteria and/orsynthesizing designs that satisfy the performance criteria.

SUMMARY

Embodiments described herein involve a method of classifying designcriteria comprising receiving design criteria for a product part. Thecriteria comprise one or both of performance and manufacturing criteria.The design criteria are sorted into different classes of one or both ofone or more objective functions and one or more constraints based onwhen they can be satisfied or optimized. Constraint violations aredetermined. A design workflow is produced to generate one or moredesigns of a part to comply with one or more of satisfying constraintsand optimizing objective functions.

Embodiments involve a system for classifying design criteria. The systemcomprises a processor and a memory storing computer program instructionswhich when executed by the processor cause the processor to performoperations. The operations comprise receiving design criteria for aproduct part. The criteria comprise one or both of performance andmanufacturing criteria. The design criteria are sorted into differentclasses of one or both of objective functions and constraints based onwhen they can be satisfied or optimized. Constraint violations aredetermined. A design workflow is set up to generate one or more designsof a part to comply with one or more of satisfying constraints andoptimizing objective functions.

A method involves receiving one or more constraints for a design of aproduct part. A feasible design space is defined based on the one ormore constraints. The design space is pruned based on a subset or all ofthe one or more constraints. The pruned design space is explored toproduce one or more feasible part designs.

A system comprises a processor, a user interface, and a memory storingcomputer program instructions which when executed by the processor causethe processor to perform operations. The operations comprise receivingone or more constraints for a design of a product part via the userinterface. A feasible design space is defined based on the one or moreconstraints. The design space is pruned based on a subset or all of theconstraints. The pruned design space is explored to produce one or morefeasible part designs.

The above summary is not intended to describe each embodiment or everyimplementation. A more complete understanding will become apparent andappreciated by referring to the following detailed description andclaims in conjunction with the accompanying drawings.

BRIEF DESCRIPTION OF THE DRAWINGS

FIGS. 1A and 1B illustrates functional load-bearing surfaces where thelatch mates with other parts in accordance with embodiments describedherein FIGS. 2 and 3 illustrate design processes for an example productpart in accordance with embodiments described herein;

FIG. 4 shows a set of design solutions on the Pareto frontier defined bymass and compliance as conflicting objectives in accordance withembodiments described herein;

FIGS. 5A-5H and 6A-6B illustrate filtering of the sensitivity field inaccordance with embodiments described herein;

FIG. 7 illustrates an example design optimization framework that may ormay not result in a feasible design in accordance with embodimentsdescribed herein;

FIG. 8 shows a design optimization framework using design space pruningand design space exploration in accordance with embodiments describedherein;

FIG. 9A illustrates a process for classifying design criteria inaccordance with embodiments described herein;

FIG. 9B shows a method for producing feasible part designs in accordancewith embodiments described herein;

FIG. 10 illustrates how solving forward problems by different FP-solversis abstracted by mapping a design (geometry and material properties) todifferent fields such as deformation, stress, and accessibility fieldsin accordance with embodiments described herein;

FIGS. 11A-11C show that design space pruning can be abstracted byintersecting the design space with feasibility halfspaces in accordancewith embodiments described herein;

FIGS. 12A-12D show that the dual properties of the FP- and IP-solvers(i.e., Sweep and Unsweep) can be leveraged to construct an exact orapproximate representation of a maximal design in accordance withembodiments described herein;

FIGS. 13A-13C shows a manufacturing setup showing six clamps used tolocate and hold a designed part and a 2-axis instrument that can move inthe plane in accordance with embodiments described herein;

FIG. 14 compares the solution designs based on a design domain to a TOproblem with and without constraining the support material volume inaccordance with embodiments described herein;

FIG. 15A shows an example Pareto tracing process in accordance withembodiments described herein;

FIG. 15B illustrates incremental steps in a Pareto tracing process inaccordance with embodiments described herein;

FIG. 16A shows a preprocesses for a latch where functional surfaces arespecified in accordance with embodiments described herein;

FIG. 16B illustrates that the Unsweep removes parts of the pre-processedinitial design that would exit the envelope for a pre-specifiedclockwise rotation of 21 degrees around the pivot in accordance withembodiments described herein;

FIG. 17 illustrates the Pareto front for solving the above problem,starting from the pruned design domain in accordance with embodimentsdescribed herein;

FIG. 18 shows problem of design for manufacturability with functionalfeatures that are specified in pre-processing in accordance withembodiments described herein;

FIG. 19 illustrates the Pareto front of accessible designs optimizedunder specified loading boundary conditions in accordance withembodiments described herein;

FIG. 20 shows the fixturing setup, raw stock, maximal manufacturabledomain (i.e., initial design for TO), and the optimized design at 40%volume fraction in accordance with embodiments described herein;

FIG. 21 shows penalizing TSF by the inaccessibility measure for T-shapedtool approaching from left or from both left and right in accordancewith embodiments described herein;

FIG. 22 illustrates optimized shapes with and without the accessibilityconstraints with one and two tool orientations in accordance withembodiments described herein;

FIG. 23 shows the Pareto fronts with one approach direction, twoapproach directions, and an unconstrained solutions in accordance withembodiments described herein;

FIG. 24 shows the optimized latches at 35% volume fraction with andwithout the accessibility constraint in accordance with embodimentsdescribed herein;

FIG. 25A shows the original TSF for compliance in accordance withembodiments described herein;

FIG. 25B illustrates the inaccessibility measure obtained from aconvolution of the design and tool at 0° and 1800 in accordance withembodiments described herein;

FIG. 25C shows the penalized TSF to incorporate accessibility and retainfunctional surfaces for the final design at a volume fraction of 35% inaccordance with embodiments described herein;

FIG. 26 illustrates an example computer architecture that can implementmethods described herein;

FIGS. 27, 28, 29, 30, 31, and 32 are formulas in accordance withembodiments described herein.

The figures are not necessarily to scale. Like numbers used in thefigures refer to like components. However, it will be understood thatthe use of a number to refer to a component in a given figure is notintended to limit the component in another figure labeled with the samenumber.

DETAILED DESCRIPTION

Mechanical design problems require reasoning about diverse, multiple,and often conflicting objectives and constraints arising fromrequirements across a product's lifecycle. The engineering designchallenge lies in traversing the trade space of these requirements tosynthesize feasible designs. This challenge has recently been amplifiedby rapid advances in manufacturing processes. Light-weight,high-performance, and multi-material composite structures with complexgeometry and material distribution can now be fabricated using variousadditive manufacturing (AM) processes. Yet, existing computer-aideddesign (CAD) systems are far behind in their representations andalgorithms to navigate the high-dimensional trade spaces that growexponentially in the number of available decisions per spatial elements.Additional functional constraints such as manufacturability, ease ofassembly, motion in presence of obstacles, and aesthetics dramaticallyincrease the trade space complexity.

Specialized domain-specific computational tools are used to generatedesigns that satisfy specific types of functional requirements. Forexample, to maximize a part's performance with as little cost ormaterial as possible, one may employ topology optimization (TO) tools.In most TO approaches, an objective function is defined over a designdomain in terms of the physical performance (e.g., strength and/orstiffness) with additional constraints on total mass or volume as wellas boundary conditions (e.g., loading and/or restraints) that oftenaccount for interfaces with other parts. TO produces valid designs withnonintuitive shapes, topologies, and material distributions that meetphysical requirements, but is rarely aware of other design criteria suchas kinematic constraints. On the other hand, to ensure collision-freemotion of a part in an assembly, one may need to examine its freeconfiguration space to guarantee collision avoidance. Similarly, forsubtractive manufacturing (SM), the machinability of a designed part ispredicated on whether the volume to be removed from a raw stock isaccessible within the cutting tool assembly's non-collidingconfigurations. For AM, one may need to consider the part's morphology,minimum feature size, and skeleton. Hybrid manufacturing (combined AMand SM) requires more complicated logical reasoning. These problems mayuse nontrivial interference analysis of shapes in relative motion thatrely on different tools of reasoning than physics-driven design toolssuch as TO. The latter often ignore motion related constraints byconsidering them out-of-scope.

Generating practical designs includes simultaneous reasoning aboutshape, motion, materials, physics, manufacturing, and assembly, amongother factors. For example, a machine part that moves relative to otherparts in a mechanical device has to avoid collisions with bothstationary and moving obstacles. These requirements are imposed askinematic constraints, expressed in terms of pointset containment ornon-interference. The same part has to sustain mechanical loads at itsjoints or points of contact with other parts. These constraints may beimposed as physical constraints, expressed in terms of (in)equalities ofmathematical functions that represent physical response (e.g., bounds ondeflection or stress). Moreover, the part has to be manufacturable withone or more AM or SM capabilities. Manufacturability constraints can beof both kinematic and physical types; for instance, accessibility in SMand post-processing of AM (e.g., support removal) are of predominantlykinematic nature, whereas achieving desired material properties in AMrequires in situ physical analysis. With few exceptions (e.g., TO for AMwith minimized support) TO algorithms are not developed withmanufacturability provisions built into their objective functions.

Different computational services, herein called solvers, are used toassist with the design process with heterogeneous constraints byproviding either analysis tools to evaluate the performance of one ormore given designs, or synthesis tools to generate one or more designsthat satisfy a given set of performance criteria analysis. These twotypes of solvers are referred to herein as forward and inverse problemsolvers (‘FP/IP-solvers’), respectively. Specifically, generative designtools are IP-solvers which solve the inverse problem by systematicallygenerating candidate designs and evaluating their performance (usingFP-solvers) to guide refinement of designs until the criteria are met.

It is unlikely that a single IP-solver is capable of simultaneousreasoning about all design criteria (e.g., objective functions andconstraints). Therefore, a typical computational design workflowrequires carefully organizing and reasoning about severalmultidisciplinary solvers to address heterogeneous design criteria.While every IP-solver reasoning about a subset of these criteria mayproduce designs that are feasible (and ideally, optimal) only withrespect to the specific criteria it considers, it provides no guaranteesabout the rest of the criteria. Different IP-solvers are thus likely togenerate designs distinct from one another, while none of themsimultaneously satisfies all criteria. Except for extremely simplecriteria, it appears impossible to combine these solutions in anyobvious way that preserves the constraints satisfied separately by eachsolver, or at least provides the best compromise. Even if such asolution exists, there may not exist any particular ordering of thesesolvers to find it, simply because each solver performs prematureoptimization with regards to the subset of criteria it cares about.

Consider the problem of designing a car hood latch. An initial designdomain with boundary conditions is provided, and the goal is to find adesign that is as stiff as possible with the least mass. Moreover, thelatch has to be free to rotate clockwise by 20 degrees around a revolutejoint, without exiting the initial design domain, so that it would notcollide with other parts that are possibly outside that envelope.Functional load-bearing surfaces where the latch mates with other partsare shown in FIGS. 1A and 1B. All feasible designs retain these surfacesas specified. The said requirements immediately suggest using twoIP-solvers that are well-positioned to deal with them (each solversatisfying a subset of them):

-   -   Let Unsweep be a solver that generates a design that remains        within a given region of space while moving according to a        prescribed motion (in this case, a clockwise rotation of 21        degrees).    -   Let PareTO be a solver that, starting from an initial design,        generates a design on the Pareto front of the two objectives        (compliance and mass), i.e., one that satisfies the stiffness        requirement and the given boundary conditions with minimal mass.

FIG. 2 illustrates a design process for an example product part inaccordance with embodiments described herein. An initial design 210 isused to create the feasible designs. Using the Unsweep and the PartTOsolvers separately, one may generate two distinct designs 220, 230, asillustrated in FIG. 2. However, there is no clear operation with whichto combine these two, in order to generate a design that satisfies bothkinematics- and physics-based constraints. For example, the intersectionof the two designs, shown in FIG. 2 generates a design 240 that does notviolate the constraints satisfied by Unsweep, because every subset ofthe unswept volume also satisfies the containment constraints. However,the constraints 242, 244 satisfied by PareTO are no longer satisfied,because the load paths are changed due to the changed topology and thecompliance target is no longer met.

In the example shown in FIG. 2, Unsweep has a property that its solutioncan be interpreted not as a single design, but as a representation ofall designs that satisfy the containment constraints. This family ofdesigns is closed under set intersection, i.e., intersecting any of thefeasible designs with another set leads to another feasible design.Similarly, PareTO can generate a family of designs that satisfycompliance requirements for different mass budgets. However, this familyis not closed under set intersection.

It is possible to obtain a feasible solution to the latch design problemby using the same set of IP-solvers, if the workflow is organizeddifferently, as shown in FIG. 3. Suppose the input design criteria 310is used to solve for the containment constraint using Unsweep. A validintermediate design 320 is generated that does not exit a given envelopewhile moving. This intermediate design 320 can be used as the initialdesign input to PareTO and optimize its shape and topology to achievethe compliance target 330 with minimal mass. This approach works becausePareTO is a material reducing solver, i.e., its solutions remainstrictly contained within the initial design. Hence any design generateddownstream will be faithful to the containment constraint that wassatisfied upstream. The same argument is not true if the order ofapplying the solvers is swapped. There is no reason to believe thatapplying Unsweep to a topologically optimized solution of PareTO willremain on the Pareto front. The fundamental differences between the twoIP-solvers should be taken into account when deciding on theirarrangement in a workflow—in this case, choosing between parallel orsequential execution and the proper order for the latter.

FIG. 4 shows a set of design solutions 410, 420, 430, 440 on the Paretofront based on an initial design 450. These design solutions 410, 420,430, 440 defined by mass and compliance as conflicting objectives.PareTO provides a clear advantage over classical TO with a singlesolution, by producing many alternatives—some of which might satisfyadditional constraints that were not accounted for in TO. Thisgenerate-and-test approach is plausible, and in fact, is a commonstrategy for dealing with heterogeneous constraints. However, it turnsout that none of the Pareto-optimized designs in this case will notremain within the envelope after a clockwise rotation of 21 degrees. Inthe example shown in FIG. 4, none of the topologically optimized designalternatives obtained from applying PareTO to the initial design 450will satisfy the containment constraint. Generate-and-test does notalways work.

In general, a good rule of thumb to organize solvers in a workflow is tocall the ones that produce the broadest families of designs earlier. Theupstream IP-solvers should generate a large number of designs, asopposed to fixing one or few designs, to provide more flexibility fordownstream solvers. The downstream solvers can be FP-solvers, testingfor new constraints, or IP-solvers, applying further optimization.However, each solver may prematurely optimize designs that may failevaluation criteria considered in downstream solvers. The “blind”process of generating and testing designs without carefully consideringproperties of the workflow- and the associated feasible design space foreach solver in the workflow—will scale poorly with increasing number ofconstraints/solvers and their complexity. A systematic approach toarranging solvers into workflows that guarantees satisfying newconstraints without violating the already satisfied constraints is muchdemanded.

Real-world design problems often involve solving multi-objectiveoptimization problems where the goal is to find the best trade-offbetween multiple competing objective functions. Classical methods suchas linear programming (e.g., the ‘simplex’ algorithm), nonlinearprogramming (e.g., the steepest descent and conjugate gradients), orNewton-Raphson are limited to single-objective optimization problems.Even for single-objective optimization, finding the global optimum isNP-hard. Numerous approaches have been developed to converge to locallyoptimal solutions in reasonable computation time to multi-objectiveproblems across different disciplines.

Unlike single-objective optimization, a total ordering for feasiblesolutions may not be possible in multi-objective optimization, i.e.,there may not exist a single “best” solution due to competingobjectives. However, feasible solutions may be partially orderedaccording to Pareto efficiency—also known as Pareto-Koopmans efficiencyor dominance. Pareto-optimal solutions are locally optimal (according toPareto-efficiency) where improving one objective comes at the expense ofat least one other objective. The collection of all Pareto-optimalsolutions is referred to as a Paretofront, which represents a curve,surface, or higher-dimensional manifold in the design space for two,three, or higher number of competing objectives, respectively. Tracing aPareto front is a challenge in multi-objective and multi-disciplinarydesign optimization.

TO has emerged as a practical class of computational methods fordesigning high-performance light-weight structures and has been appliedin numerous areas such as designing automobiles components, aircraftcomponents, spacecraft modules, cast parts, compliant mechanisms, andmany other products. Numerous approaches such as density-based,levelset-based, and evolutionary methods for TO have been developed.

TO typically focuses on optimizing designs for performance (e.g.,physical response to loads during operation) but less on other factorssuch as manufacturability. Apart from traditional processes such asmachining and molding, more recent technologies such as AM haveintroduced the ability to fabricate complex topologically optimizeddesigns while presenting new manufacturing challenges Processlimitations may be considered during the design/optimization stage asmuch as possible to avoid repeated prototyping and iterations until theoptimized designs are manufacturable. Specifically, applying correctionsto the geometry or topology of a solution after TO to make itmanufacturable may sacrifice the achieved optimality.

One solution is to impose design rules obtained from domain expertiseand experience. These rules relate specific combinations of shape,materials, and process to impose simplified constraints that can bebuilt into the TO framework to restrict the feasible design space. Forexample, when designing for AM via fused deposition modeling usingpolymers, one should require that all facets oriented at an anglegreater than 45 degrees with respect to the build direction be supportedwith additional scaffolding material. When designing for casting andinjection molding processes, one should ensure that the part hasfeatures of almost uniform thickness and no entrapped holes are present,so that the mold can be removed and the molten material cools downuniformly throughout the part. When designing for wire- orlaser-cutting, one should ensure that the final design has a uniformcross-section, i.e., is 2.5D along the cutting direction. Theseconstraints can be imposed during TO through filtering of thetopological sensitivity field (TSF) as illustrated in FIGS. 5A-5H.

FIG. 5A illustrates the initial design domain having boundary conditionsas shown. FIG. 5B shows a resulting optimized design withoutmanufacturing restraints. According to embodiments described herein, thedesign shown in FIG. 5B has no TSF filtering. FIGS. 5C and 5D show anoptimized design with wire and/or laser-cut constraints with through-cutfilter 520 having unfiltered TSF 510 as shown. This involves designing apart to have a 2.5D geometry, i.e., have a fixed 2D section along thewire or laser beam direction 530. FIGS. 5E and 5F illustrate anoptimized design with casting constraints, that involves designing apart to have a monotonically reducing 2D section along the drawdirection. FIGS. 5G and 5H show an optimized design with casting(similar to FIG. 5F) and additional constraints that involve retainingsome functional surfaces. The different filters for FIGS. 5C-5H areshown as curves along the arrow shown on a cross-section of the part.Each filter modifies the unfiltered TSF such that its level setsatisfies the design constraints, producing different optimized designs.

Another AM consideration is the manufacturing resolution, which can bedirectly incorporated into the TO algorithm as a minimum feature sizeconstraint through either local gradient constraints or TSF filtering asshown in FIGS. 6A and 6B. FIG. 6A illustrates a first design 610. TSFfiltering is applied 620 that constrains the minimum feature size of thedesign. A second design 630 is generated with the minimum feature sizefiltering applied. FIG. 6B also illustrates two curves that representthe cross-sections of the unfiltered TSF 650 (with a smallerpeak-to-peak distance of 5 mm) and a filtered TSF 640 (with a largerpeak-to-peak distance of 10 mm), where the latter would produce a designwith a larger minimum feature size as the level set of the filtered TSF.

It is also possible to reduce the amount of support structure that maybe needed in an AM by either finding a good build orientation or TSFfiltering. Build orientation optimization often involves solving amulti-objective problem taking into account other factors such assurface quality, build time, or manufacturing error. TSF filtering, onthe other hand, can be achieved by penalizing overhang surfaces,penalizing undercut surfaces, and/or augmenting new TSFs.

Solving heterogeneous (e.g., kinematic, physical, and manufacturing)constraints with multidisciplinary solvers is difficult. Each solvermight make decisions with care for its target subset of constraintswhile potentially violating the rest of the constraints. FIG. 7illustrates an example design optimization framework that may or may notresult in a feasible design. Arranging multiple solvers sequentially maywork in some orders and fail in others, depending on what propertiesthey preserve. The first example 700 shows the solvers Unsweep 710 andPareTO 720 called in series with the Unsweep solver being called firstas described previously in conjunction with FIG. 3. According to variousconfigurations, one or more other solvers 730 may be used to generate afinal feasible design. The second example 735 again shows the solversUnsweep 750 and PareTO 740. In this case, the PareTO solver 740 iscalled first. One or more additional solvers 760 may be called dependingon the given constraints. In this example, there are no feasibledesigns.

The third example, 765 shows the Unsweep solver 770 and the PareTOsolver 780 called in parallel. The results of the Unsweep solver 770 andthe PareTO solver 780 are then combined 785 as described in conjunctionwith FIG. 2. One or more additional solvers 790 may be called dependingon the given constraints. The order of this example may result infeasible designs. Parallel composition requires combining solutions inways that do not always preserve the properties either. A two-phaseapproach that includes design space pruning and design space explorationmay be used to create feasible designs. The former invokes IP-solvers inan arbitrary order to cut out the infeasible design subspace withoutpremature optimization. The latter navigates the pruned design space byfixed-point iterations over FP-solvers.

In addition or as an alternative to the above examples of performanceand manufacturing requirements, there may be other design criteria thatinvolve spatial reasoning about the interactions of moving (translatingand rotating) shapes such as collision avoidance, packaging, robotmotion planning, and accessibility analysis. These requirements cannotbe easily enforced by design rules, TSF filtering, or other techniquescommonly used in TO. Rather, they are often expressed as setconstraints, i.e., statements in the language of sets (e.g., in terms ofaffine transformations, Boolean operations, and containment) rather thanthe language of real-valued functions used for (in)equality constraintsin TO. A broad class of inverse problems in practical design andmanufacturing reduce to solving set constraints formulated in theconfiguration space of rigid motions.

Although the problems with set constraints are common, they are notmainstream in design/optimization workflows due to non-smoothness andcomputational intensity. There are instances of TO frameworks that dealwith motion-related problems in an ad hoc manner; for instance, inmodeling collision and contact when designing compliant mechanismsand/or parts made of hyperelastic materials that undergo largedeformations part. It may not be immediately obvious how set constraintscan be incorporated in a systematic fashion into the design/optimizationprocess without incurring prohibitive computation costs of spatialanalysis at every iteration. Some examples present a differentclassification of constraints that enables design space pruning andexploration, in which set constraints are also restated in terms of(in)equalities of functions.

According to embodiments described herein, solvers that reduce materialfrom a bounded domain in 2D and 3D space are used to generate designs.When all objectives and constraints cannot be handled by a singlesolver, current practice relies on a case-by-case domain-specificanalysis to properly construct workflows FIG. 7 depicts examples ofcombining Unsweep and PareTO. An arrangement that works in one case maynot work in another.

Solvers described herein may be classified into two fundamental types:

1. Design space pruning solvers restrict the feasible design space bypruning the subspaces that violate one or more design criteria. They arepermutative, meaning that they can be called at the beginning of thedesign workflow in an arbitrary order. By directly operating on thedesign subspaces as first-class entities, they postpone optimization todownstream solvers; and

2. Design space exploration solvers simultaneously explore the (pruned)design subspaces for optimized solutions.

FIG. 8 shows a design optimization framework using design space pruningand design space exploration in accordance with embodiments describedherein. Initial design constraints 810 are input into the framework. Oneor more design space pruning solvers 820, 830, 840 are called in seriesin an arbitrary and each output an optimized design based 825, 835, 845on the previous solvers optimized design. Once the pruning stage iscomplete, one or more design space exploration solvers 850, 860, 870,880 are used to explore the PareTO subspace to create one or more finaloptimized designs.

FIG. 9A illustrates a process for classifying design criteria inaccordance with embodiments described herein. Design criteria for aproduct part are received 910. According to various implementations, thedesign criteria are received via a user interface, for example. Thedesign criteria may comprise one or both of performance andmanufacturing criteria. The design criteria are sorted 920 intodifferent classes of one or both of one or more objective functions andone or more constraints based on when they can be satisfied oroptimized. Constraint types comprise one or more of set constraints,equality constraints, and inequality constraints. According toembodiments described herein, constraint scopes comprise one or more ofglobal, local, and strictly local. Constraint violations are determined930. A design workflow is produced 940 to generate one or more designsof a part to comply with one or more of satisfying constraints andoptimizing objective functions.

Various configurations involve computing one or more performance fieldsof the one or more designs. The one or more objective functions areevaluated functions based on the performance fields. The one or moredesigns are ordered based on the evaluated objective functions.According to various embodiments, one or more constraint violations arecomputed based on the performance fields and it is determined whetherthe one or more designs are feasible based on the computed constraintviolations.

FIG. 9B shows a method for producing feasible part designs in accordancewith embodiments described herein. One or more constraints for a designof a product part are received 950. According to variousimplementations, the one or more constraints are received via a userinterface. A feasible design space is designed 960 based on the one ormore constraints. The design space is pruned 970 based on a subset orall of the one or more constraints. The pruned design space is explored980 to produce one or more feasible part designs.

Various configurations involve evaluating point membership tests basedon one or more feasibility predicates defined by pointwise constraintsand computing maximal designs to represent the feasible design spacewith respect to the pointwise constraints. According to variousembodiments, a measure of global changes of each of a plurality ofconstraint violations due to local variations in the one or more designsis computed. A topological sensitivity field (TSF) may be determined bycombining the measures of global changes of each of the plurality ofconstraint violations, the TSF configured to guide design spaceexploration to produce the feasible designs. According to variousimplementations, defining the feasible design space comprises using oneor more of a density-based, a levelset shape, and a topologyoptimization approach. The one or more designs may be optimized usingaugmented Lagrangians.

FP-solvers (e.g., FEA and manufacturability analysis) map an instance ofthe design space (i.e., a “design”) to an instance of the performancespace (i.e., a “field-tuple”). Predicates are defined to decide whetherthe design is satisfactory with respect to constraints in terms ofperformance variables. An IP-solver modifies the design until allconstraints are satisfied.

FIG. 10 illustrates how solving forward problems by different FP-solversis abstracted by mapping a design (geometry and material properties) todifferent fields such as deformation, stress, and accessibility fieldsin accordance with embodiments described herein. Each field is evaluatedagainst the constraints, whose satisfaction are captured by the binarypredicates. Inverse problem solving by IP-solvers involves updating thedesign until all constraints are satisfied.

The ‘performance’ of a given design Ω∈D is an n-tuple

(Ω):=(

₁(Ω), . . . ,

_(n)(Ω)). Think of the performance space as a product space P:=(F₁× . .. ×F_(n)), where each F_(i) is a class of fields, i.e., each

_(i)∈F₁ is an integrable field

_(i)(Ω): Ω₀→

_(i) over the design domain Ω₀ whose value at a given “query point” x∈Ω₀is denoted by

_(i)(x; Ω):=(

_(i)(Ω))(x)∈

_(i). Examples of such fields are:

-   -   binary-valued fields (        _(i):={0,1}), used to describe indicator functions of regions of        interest within the design domain such as non-manufacturable        features or regions requiring design correction.    -   integer-valued scalar fields (        _(i):=Z), used to characterize local topological properties of        3D printed parts or to classify atomic units of manufacturing in        hybrid (combined AM and SM) manufacturing.    -   real-valued scalar fields (        _(i):=R), 3D vector fields (        _(i):=R³), and higher-rank tensor fields used to represent        distributed physical quantities such as displacement, velocity,        stress, strain, and so on, or manufacturability measures.

Forward problem solvers (FP-solvers) map a given design instance to oneor more performance fields, hence can be viewed as implementations ofone or more maps

_(i):D→P_(i). The entire forward problem, solved by one or moreFP-solvers, can be viewed as a single map from design space toperformance space

:D→P, which has a unique outcome for a given design.

For example, consider a finite elements analysis (FEA) FP-solver thatcomputes (discretized forms of) a displacement field u_(Ω):=

₁(Ω) for small deformations of a given design Ω∈D due to boundaryconditions such as restraints and external forces. It also may beimportant to compute the stress field σ_(Ω):=

₂(Ω), which depends locally on displacement and material properties(e.g., the linear elasticity law). The vector/tensor values of solutionfields probed at a query point x∈Ω₀ are denoted by u_(Ω)(x)=

₁(x; Ω) and σ

(x)=

₂(x; Ω). In this case, both functions are zero outside the design, i.e.,

_(1,2) (x; Ω)=0 if x∈(Ω₀−Ω). FEA solves the weak form of the governingdifferential equation, discretized into a linear system (e.g., using hatfunctions) [K_(Ω)][u_(Ω)]=[f], where the stiffness matrix [K_(Ω)] andexternal load vector [f] depend on the design shape and materialproperties, as well as boundary conditions. The equations are solved toobtain the discrete form of the displacement field [u_(Ω)] from whichthe discrete form of the stress field [an] is computed by linearoperations.

Another example is accessibility analysis for machining (e.g., millingor turning). For instance, consider an FP-solver for 3-axis millingsimulation, which computes (discretized forms of) a volumetric measureof inaccessibility as a field μ_(Ω)(x):=

₃(Ω) for a given design Ω∈D and machine tool parameters. This measure ata query point x∈Ω₀, denoted by μ_(Ω)(x)=

₃(x; Ω), returns the penetration volume of the moving tool assemblyT=(H∪C), including the holder H and cutter C, into the stationaryobstacles O_(Ω)=(Ω∪F), including target form Ω and fixtures F. Thesolver computes the discrete form of this field (e.g., sampled at pointclouds or voxels) [μ_(Ω)] as a convolution between the discrete forms ofthe indicator functions of stationary solids [1_(O) _(Ω) ] and movingsolids [1_(T)] using a fast Fourier transform (FFT). The maximal set ofaccessible configurations—in this case, pure translations in 3D for afixed orientation—is then obtained as the null set M_(Ω)=μ_(Ω) ⁻¹(0),i.e., the translations that do not lead to undesirable collisions. Themaximal removable volume is obtained by sweeping the cutter with themaximal motion, i.e., R_(Ω):=sweep(M_(Ω), C). Its indicator function1_(R) _(Ω) (x)=

₄(x; Ω) can be viewed as a predicate for accessibility, i.e., it returns1 (resp. 0) if the query point x∈Ω₀ is (resp. is not) accessible. Thediscrete form of this binary field [1_(R) _(Ω) ] can also be obtained bythresholding FFT-based convolution of discrete forms [1_(M) _(Ω) ] and[1_(C)].

The computations performed by the above two FP-solvers (FEA andaccessibility analysis) are abstracted by

:=(

₁,

₂,

₃,

₄) that maps a given design to a “field tuple” that represents analysisresults. FIG. 10 illustrates one instance of each such field for atopologically optimized bracket.

Inverse problem solvers (IP-solvers), on the other hand, find one ormore designs that satisfy a given collection of functional requirements.Most IP-solvers employ an iterative process to:

1. generate one or more valid candidate design(s);

2. perform analysis on the candidate design(s) to compute theperformance(s) of interest (using one or more FP-solver(s));

3. evaluate the performance(s) against given functional requirements;and,

4. if the requirements are not met, decide on the next generation ofcandidate design(s) based on the current evaluation and update rules.

The process is repeated until the requirements are met. The evaluationprocess (item 3) can be conceptualized a finite number of predicatesdefined over the performance space as maps c_(i):D→{0,1} for i=1, 2, . .. , n. Each predicate's outcome c_(i)(Ω)∈{0,1} is determined by means ofa constraint imposed on the performance field

_(i)(Ω)∈F_(i) simulated by an FP-solver as shown in (1). Without loss ofgenerality, it can be assumed t every requirement depends on oneperformance field only. If more than one field is used in calculating apredicate, those fields can be tupled into another field. If more thanone requirement is computed on one field, it can be thought of as twocopies of the same field.

$\begin{matrix}{{c_{i}(\Omega)}:=\left\{ \begin{matrix}{{1\mspace{14mu}{if}\mspace{14mu}{the}\mspace{14mu} i^{th}\mspace{14mu}{constraint}\mspace{14mu}{is}\mspace{14mu}{satified}},} \\{0\mspace{14mu}{{otherwise}.}}\end{matrix} \right.} & (1)\end{matrix}$

These can be (in)equality constraints, which are common in physics-baseddesign formulations as in TO, and set constraints, which are ubiquitousin design under kinematics-based constraints such as packaging,assembly, and accessibility for manufacturing.

According to embodiments described herein, the performance criteriaevaluation can be thought of as a map c:=(c₁, c₂, . . . ,c_(n)):D→{0,1}^(n) i.e., c(Ω) is a binary string whose bits indicatewhether a given design satisfies each of the criteria. A design Ω∈D iscalled ‘feasible’ if it simultaneously satisfies all criteria, i.e.,c(Ω)=(1, 1, . . . , 1). The feasible design subspace D*⊆D is the subsetof all feasible designs, defined as in (2).

D*:={Ω∈D|c(Ω)=1^(n) }=:c ⁻¹(1).  (2)

Here, (Ω)⁻¹ denotes inversion of a mathematical relation. Given ƒ:X→Y,ƒ⁻¹(y)={x∈X|ƒ(x)=y}. Note that in general, ƒ⁻¹(y)⊆X is a set, i.e.,ƒ⁻¹:Y→

(X). It is a singleton set (i.e., has one element) if the function isbijective.

Unlike forward problems, inverse problems have non-unique, ofteninfinitely many, solutions (i.e., |D*|>1). The feasible design space canalso be defined as the intersection of the feasibility halfspacesH_(i):=c_(i) ⁻¹(1), each implicitly describing one of the designsubspaces that satisfy one criterion at-a-time as shown in (3).

D*=∩ _(1≤i≤n) c _(i) ⁻¹(1)=D−∪ _(1≤i≤n) c _(i) ⁻¹(0).  (3)

The idea of design space pruning is to progressively cut out portions ofthe design space that violate any one of the criteria. Theoretically,pruning can be done in an arbitrary order-noting that intersections orunions of the halfspace in (3) are permutative. Computationally,however, it is only possible if the design subspaces H_(i)=c_(i) ⁻¹(1)can be manipulated by algorithms as first-class entities. The goal is tounderstand under what conditions such manipulations are possible and howan entire design subspace can be represented.

The predicates introduced earlier are implemented in practice by testingwhether a candidate design's performance satisfies an (in)equalityconstraint. Note that every equality constraint g(⋅)=0 can berepresented by two inequality constraints g(⋅)≤0 and −g(⋅)≤0. Suchconstraints can be classified into three types; namely, global, local,and strictly local (in)equality constraints.

It is common to have design criteria specified in terms of globalconstraints g_(i)(Ω)≤0, i.e., by defining a predicate of the generalform as shown in (4).

$\begin{matrix}{{c_{i}(\Omega)}:=\left\{ \begin{matrix}{{{1\mspace{14mu}{if}\mspace{14mu}{g_{i}(\Omega)}} \leq 0},} \\{0\mspace{14mu}{{otherwise}.}}\end{matrix} \right.} & (4)\end{matrix}$

Here, g_(i):D→R is a function of the entire shape of the design Ω∈D,potentially in addition to fixed external factors such as boundaryconditions, manufacturing process parameters, packaging envelope,operating conditions (e.g., motion in assembly), etc. The constraint isoften evaluated in terms of a global property of an entire performancefield

_(i)(Ω)∈F_(i), e.g., as an upper/lower-bound on its maximum/minimum orits integral properties such as p-norms. This is denoted by g_(i)(Ω)=g_(i)(

_(i)(Ω)) where g _(i):D→R. For example, one can constrain the maximaldisplacement or maximal stress of a solid under external loads (FIG.3.1) by using the following constraints in (4):

$\begin{matrix}{{{g_{1}(\Omega)} = {{{\overset{\_}{g}}_{1}\left( u_{\Omega} \right)}:={{\max\limits_{x \in \Omega}{{u_{\Omega}(x)}}} -_{UB}}}},} & (5) \\{{{g_{2}(\Omega)} = {{{\overset{\_}{g}}_{2}\left( \sigma_{\Omega} \right)}:={{\max\limits_{x \in \Omega}{{\sigma_{\Omega}(x)}}} - \sigma_{UB}}}},} & (6)\end{matrix}$

Here, _(UB), σ_(UB)>0 are constant upper-bounds on the magnitude of thedisplacement vector and stress tensor, captured by the constraintsg₁(Ω)≤0 and g₁(Ω)≤0, respectively. One can in general use the p-norm ofthe fields for finite (but large) p≥1, noting that maximum is a specialcase as p→∞. This is especially useful to smooth out the possiblesingularities (e.g., infinite stress, due to stress concentrations).

According to embodiments described herein, (7) is another example toconstrain a design to be manufacturable via machining, using theaccessibility analysis mentioned earlier.

g ₃(Ω)= g ₃(1_(R) _(Ω) ):=∫_(Ω) ₀ _(−Ω)¬1_(R) _(Ω) (x)dv[x]−V_(UB),  (7)

Here, ¬1_(R) _(Ω) (x)=1−1_(R) _(Ω) (x) is a negation. This constraintrestricts the total volume of inaccessible regions R_(Ω)⊆Ω^(c), obtainedas the 1-norm of their indicator function, by an upper-bound V_(UB)>0.

According to embodiments described herein, it is sometimes possible todefine a predicate in terms of local constraints evaluated at a specificpoint in the design domain; for instance, using one or both of (8) and(9).

$\begin{matrix}{{c_{i}(\Omega)}:=\left\{ \begin{matrix}{{1\mspace{14mu}{if}\mspace{14mu}{\forall{x \in {\Omega_{0}:{{g_{i}\left( {x;\Omega} \right)} \leq 0}}}}},} \\{0\mspace{14mu}{{otherwise}.}}\end{matrix} \right.} & (8) \\{{c_{i}(\Omega)}:=\left\{ \begin{matrix}{{1\mspace{14mu}{if}\mspace{14mu}{\exists{x \in {\Omega_{0}:{{g_{i}\left( {x;\Omega} \right)} \leq 0}}}}},} \\{0\mspace{14mu}{{otherwise}.}}\end{matrix} \right.} & (9)\end{matrix}$

Note that the two alternative forms are different by the “for all” and“there exists” quantifiers, which may lead to different globalimplications. Unlike the case with (4), here g_(i):(Ω₀×D)→R is field fora fixed design Ω∈D, i.e., is also a function of the query point x∈Ω₀. Inturn, the constraint is evaluated based on the probed value of theperformance field

_(i)(x; Ω)∈

_(i) (generally, a tensor) at the query point. This dependency isdenoted by g_(i)(x; Ω)=g _(i)(

_(i)(x; Ω)) where g _(i):

_(i)→R. For example, the global displacement and stress bounds mentionedearlier can be imposed locally as shown in (10) and (11).

g ₁(x;Ω)= g ₁(u _(Ω)(x)):=∥u _(Ω)(x)∥_(UB),  (10)

g ₂(x;Ω)= g ₂(u _(Ω)(x)):=∥σ_(Ω)(x)∥−σ_(UB),  (11)

It is easy to verify that using (10) and (11) with (8) is equivalent tousing (5) and (6) with (4) in this example. In general, localconstraints g_(i)(x; Ω)≤0 used with “for all” or “there exists”quantifiers in (8) and (9) can be equivalently expressed as globalconstraints (stated independently of x∈Ω₀) via min/max, respectively:

$\begin{matrix}{\left. \left\lbrack {\forall{x \in {\Omega_{0}:{{g_{i}\left( {x;\Omega} \right)} \leq 0}}}} \right\rbrack\rightleftharpoons{\max\limits_{x \in \Omega_{0}}{g_{i}\left( {x;\Omega} \right)}} \leq 0 \right.,} & (12) \\\left. \left\lbrack {\exists{x \in {\Omega_{0}:{{g_{i}\left( {x;\Omega} \right)} \leq 0}}}} \right\rbrack\rightleftharpoons{\min\limits_{x \in \Omega_{0}}{g_{i}\left( {x;\Omega} \right)}} \leq 0. \right. & (13)\end{matrix}$

As another example, consider the accessibility analysis discussedearlier. Instead of constraining the total volume of inaccessibleregions via the global constraint of (7), the inaccessibility measurecan be locally constrained as shown in (14).

g ₄(x;Ω)= g ₄(μ_(Ω)(x)):=(1_(O) _(Ω) ,*1_(−T))(x)−μ₀,  (14)

Here, O_(Ω)=(Ω∪F) and T=(S∪C) are the stationary and moving solids,respectively. μ₀>0 is a small constant to provide allowance fornumerical errors. The convolution operator * is defined as:

(1_(O) _(Ω) *1_(−T))(x)=∫_(Ω) ₀ 1_(Ω)(x′)1_(−T)(x−x′)dv[x′],  (15)

Here, 1_(−T)(x)=1_(T)(−x) is a reflection with respect to the origin,hence 1_(−T)(x−x′)=1_(T)(x′−x) is the indicator function of the movingobject (i.e., tool assembly), translated to the query point x∈Ω₀. Theintegral is nonzero at integration points x′∈Ω₀ that belongs to theinterference of the translated object with the stationary obstacles.

It is not always possible to convert global constraints to localconstraints or vice versa, without defining new performance variables,e.g., in terms of the norms of existing performance fields. When

_(i)(x; Ω) is decidable independently of Ω∈D, the above two constraintslead to maximal/minimal feasible designs (in set-theoretic terms).

A special case of (8) occurs if the predicate is decidable withoutapriori knowledge of the design itself. In other words, the constraintscan be evaluated purely from a knowledge of the query point's positionx∈Ω and external factors, if any (e.g., a known rigid body motionapplied to the entire design). The predicate's result can be obtainedwithout knowing of the overall shape of the design. This is the case if

_(i)(x; Ω)=

_(i)*(x), i.e., the forward problem's solution can be evaluatedpointwise—emphasized by the overline notation. The correspondingconstraint g_(i)*(x):=g _(i)(

_(i)*(x))≤0 is hereafter called a strictly local (i.e., pointwise)constraint. The predicates in (8) or (9) in this case are dependent on Ωonly due to the logical quantifiers for the pointwise testing as shownin (16) and (17).

$\begin{matrix}{{c_{i}(\Omega)}:=\left\{ \begin{matrix}{{1\mspace{14mu}{if}\mspace{14mu}{\forall{x \in {\Omega:{{g_{i}^{*}\left( {x;\Omega} \right)} \leq 0}}}}},} \\{0\mspace{14mu}{{otherwise}.}}\end{matrix} \right.} & (16) \\{{c_{i}(\Omega)}:=\left\{ \begin{matrix}{{1\mspace{14mu}{if}\mspace{14mu}{\exists{x \in {\Omega:{{g_{i}^{*}\left( {x;\Omega} \right)} \leq 0}}}}},} \\{0\mspace{14mu}{{otherwise}.}}\end{matrix} \right.} & (17)\end{matrix}$

Hence, one can define pointwise predicates in this case by (22)

$\begin{matrix}{{c_{i}^{*}(\Omega)}:=\left\{ \begin{matrix}{{{1\mspace{14mu}{if}\mspace{14mu}{g_{i}^{*}(x)}} \leq 0},} \\{0\mspace{14mu}{{otherwise}.}}\end{matrix} \right.} & (18)\end{matrix}$

According to embodiments described herein, the pointwise predicatedefines a point membership classification (PMC) that implicitlydetermines the entire feasibility halfspace H_(i)=c_(i) ⁻¹(1) using its“representative” maximal/minimal feasible design Ω_(i)*:=c_(i) ⁻¹(1)using (16) or (17), respectively.

The physics-based constraints exemplified earlier might reduce topointwise constraints in rare examples—e.g., the stress tensor σ*(x) forhydrostatic pressure in a liquid container at rest depends on the querypoint's depth from the surface, but not on the container's designedshape. Nevertheless, pointwise constraints are ubiquitous inkinematics-based constraints that are central to applications rangingfrom assembly and packaging to manufacturing.

Many kinematics-based design criteria lead to constraints expressed inthe algebra of sets. A common form of set constraints is in terms ofcontainment: Γ(Ω)⊆E (for a fixed envelope E⊆R^(d)). The exact sameconstraint can be written in terms of non-interference: (Γ(Ω)∩O)=∅ (fora fixed obstacle O⊆R^(d)), where Γ:D→

R^(d)) is a set transformation and E:=O^(c) (i.e., complement of O).

At a first glance, these constraints appear to have a completelydifferent form than the inequality constraints described previously inthe algebra of fields. Here, it is shown that set constraints may alwaysbe reformulated as (global or local) inequality constraints by virtue ofdescribing sets with their indicator functions. However, converting themto strictly local (i.e., pointwise) constraints is possible undercertain conditions.

According to embodiments described herein, for every solid Ω∈D, itsindicator (i.e., characteristic) function 1_(Ω):Ω₀→{0,1} is defined by(19).

$\begin{matrix}{{1_{\Omega}(x)}:=\left\{ {\begin{matrix}{{{1\mspace{14mu}{if}\mspace{14mu} x} \in \Omega},} \\{0\mspace{14mu}{{otherwise}.}}\end{matrix},{i.e.},{\Omega = {1_{\Omega}^{- 1}{(1).}}}} \right.} & (19)\end{matrix}$

Hence, every containment constraint is restated as an inequalityconstraint of the form used in (8) as shown in (20).

Γ(Ω)⊆E

[∀x∈R ^(d):1_(Γ(Ω))(x)≤1_(E)(x)],  (20)

i.e., using the standard form 1_(Γ(Ω))(x)−1_(E)(x)≤0.

The above inequality constraint of (20) can also be rewritten as aglobal constraint of the form used in (4) by upper-bounding the maximumas shown in (21).

$\begin{matrix}{{{\Gamma(\Omega)} \subseteq \left. E\rightleftharpoons{\max\limits_{x \in R^{d}}\left\lbrack {{1_{\Gamma{(\Omega)}}(x)} - {1_{E}(x)}} \right\rbrack} \leq 0 \right.},} & (21)\end{matrix}$

Notice that no assumption is made regarding the properties of thepointsets Γ(Ω), E⊆R^(d). In most practical scenarios, both are solidswithin a bounded domain, taken as the design domain Ω₀∈D. In such cases,the properties of the mapping and envelope can be exploited to computethe maximum in (21) efficiently.

Moreover, if the set constraint (Γ(Ω)∩E^(c))=∅ is in terms ofregularized intersection, it can be rewritten as a global (in)equalityconstraint in terms of the volume vol(Γ(Ω)∩E^(c)]=0 (or ≤0, forconsistency), which, in turn, is an inner product of indicatorfunctions:

$\begin{matrix}{{{vol}\left\lbrack {{\Gamma(\Omega)}\bigcap E^{c}} \right\rbrack} = \left\langle {1_{\Gamma{(\Omega)}},{⫬ 1_{E}}} \right\rangle} & (22) \\{= {\int_{R^{d}}{⫬ {1_{\Gamma{(\Omega)}}(x)1_{E}d{{v\lbrack x\rbrack}.}}}}} & (23)\end{matrix}$

Further, if Γ(Ω) is a rigid transformation of Ω, the inner product turnsinto a convolution of 1_(Ω) and ¬1_(E) over the configuration space ofmotions.

Let us consider the inaccessibility analysis one more time. Formanufacturing with precision requirements (e.g., for assembly/fit), theinaccessible regions Γ(Ω):=(Ω₀−Ω)−R_(Ω) can be restricted to becompletely contained within a tolerance zone E⊆R^(d), hence formulatingthe problem as a set constraint Γ(Ω)⊆E. It has been shown that:

$\begin{matrix}{{\Gamma(\Omega)} = {\left( {\Omega_{0} - \Omega} \right) - {{sweep}\mspace{14mu}\left( {M_{\Omega},C} \right)}}} & (24) \\{= {\left( {\Omega_{0} - \Omega} \right) - \left( {M_{\Omega} \oplus C} \right)}} & (25) \\{= {\left( {\Omega_{0} - \Omega} \right) - \left( {\left( {\Omega \oplus \left( {- T} \right)} \right)^{c} \oplus C} \right)}} & (26) \\{{= {\left( {\Omega_{0} - \Omega} \right) - \left( {\left( {\Omega^{c} \ominus \left( {- T} \right)} \right) \oplus C} \right)}},} & (27)\end{matrix}$

Here, the operators (⋅)^(c),∪,∩,−,⊕,⊖ are all regularized to ensuretheir algebraic closure within the design space D and where T=(H∪C)represents the tool assembly, including the holder H and the cutter C.It is understood that Minkowski sums in (26) can be obtained from the0-superlevel set of the convolution fields of the participating sets. Ingeneral, Y=(X₁⊕X₂)

1_(Y)=sign(1_(X) ₁ *1_(X) ₂ ). More precisely:

1_(Γ(Ω))=¬_(Ω)−sign(¬sign(1_(Ω)*1_(−T))*1_(C)),  (28)

Here, * is the convolution operator in R^(d) and sign(x)=1 (resp. 0) ifx>0 (resp. x≤0) is the sign function. The latter is may be needed toconvert the real-valued convolutions to binary-valued indicatorfunctions.

In summary, set constraints, as defined here in terms of containment ornon-interference, may always be restated as global or local inequalityconstraints. Here, the conditions under which set constraints can beconverted to strictly local (i.e., pointwise) constraints is shown, toenable design space pruning.

Depending on the properties of Γ:D→

(R^(d)), the inequality constraint 1_(Γ(Ω))(x)−1_(E)(x)≤0, used inglobal or local forms of (20) and (21), respectively, may or may not berestated as a strictly local (i.e., pointwise) constraint. The goal inthis here is to articulate the conditions under which this is possible.

To enable pointwise formulation, the dependency of the PMC for Γ(Ω) on Ωis eliminated so that 1_(Γ(Ω))(x) on the left-hand side of theinequality constraint in (20) depends only on the query point x∈Ω₀ andthe fixed envelope E⊆R^(d). This can be done if the set transformation Γis itself a pointwise transformation, meaning that it can be defined byextending a transformation of 3D points γ: (R^(d) or Ω₀)→

(R^(d)) to a transformation of 3D pointsets Γ:D→

(R^(d)) by simply applying the former to every point of the pointset andunifying the results:

Γ(Ω):=∪_(xεΩ)γ(x)={x′|X′∈γ(x),x∈Ω}.  (29)

Note that γ(x) is itself a pointset, not a point, to capture the mostgeneral case. For example, it can be a curve segment or surface patchrepresenting the 1D or 2D trajectory of a point under a given one- ortwo-parametric motion, respectively.

The above refactoring is possible for many applications. For example, ifthe design has to move (when deployed in assembly) according to a knownmotion set M⊆SE(3) without exiting a safe region of space E⊆R^(d) thatcontains no obstacles, the above constraint can be used withΓ(Ω):=sweep(M, Ω), where

sweep(M,Ω)=∪_(τΣM) τΩ={τxτ∈M,x∈Ω},  (30)

This is the sweep of the designed part as it travels by the given motion(known a priori). In this case the sweep indicator function 1_(Γ(Ω)) canbe directly obtained as follows:

1_(Γ(Ω))(x)=1iff∃τ∈M:x∈τΩ, i.e., τ⁻¹ x∈Ω,  (31)

According to embodiments described herein, a PMC test for Γ(Ω) can beobtained by applying the inverse motion M⁻¹={τ⁻¹|τ∈M} to the query pointand checking if it intersects the design. The inequality constraint in(21) can thus be computed rapidly by sampling query points in the designdomain and testing intersections of their inverse trajectories with thedesign. The computation for one point does not necessarily need one tohave explicit knowledge of the results for other points. Other thanperfect parallelization (e.g., on GPU), this property enables pruningthe design space—leading to development of the Unsweep solver discussedearlier—before optimizing the design for other criteria.

In general, for pointwise transformations as in (29), the globalconstraint Γ(Ω)⊆E can be restated as shown in (32).

Γ(Ω)⊆E

[∀x∈Ω:γ(x)⊆E],  (32)

i.e., Ω remains within E after a Γ-transform iff all points inside itremain within E after a γ-transform. Note also that inequalityconstraint in (24) can now be rewritten in a pointwise fashion every x∈Ωas shown in (33).

∀x∈Ω:∀x′∈R ^(d):1_(γ(x))(x′)≤1_(E)(x′),  (33)

As with other inequality constraints, not every global or local setconstraint can be converted to a pointwise set constraint. For example,the toleranced accessibility constraints Γ(Ω)⊆E forΓ(Ω):=(Ω₀−Ω)−sweep(M_(Ω), C) in (24) cannot be evaluated pointwise,because the maximal collision-free motion M_(Ω)=(Ω⊕((−T))^(c) depends onthe global shape of Ω, unlike the case with the fixed motion in theearlier example with Unsweep.

According to embodiments described herein, a fairly general formulationof a design problem subject to n≥1 heterogeneous (e.g., kinematics- andphysics-based) constraints is shown.

Without loss of generality, let n_(C)=n_(G)+n_(L)+n_(p) where 0≤n_(G),n_(L), n_(P)≤n_(C) are the number of global, local, and strictly local(i.e., pointwise) constraints, respectively. All constraints, includingset constraints, are expressed as inequality constraints for uniformity.The design problem is to identify the feasible design space D*=c⁻¹(1),stated as a constraint satisfaction problem:

-   -   Find D*⊆D, such that for all x∈Ω∈D*:

g _(i)*(x)≤0, for 0<i≤n _(P),

g _(i)(Ω)≤0, for n _(P) <i≤n _(P) +n _(G),

g _(i)(x;Ω)≤0, for n _(P) +n _(G) <i≤n _(C),  (34)

It is assumed that none of the g_(i)(x; Ω)≤0 can be simplified into oneof g_(i)(x)≤0 or g_(i)(Ω)≤0 forms. Hereafter, the P-, G-, andL-subscripts are used for various notions related to pointwise, global,and local constraints, respectively; for instance,D*=D_(P)*∩D_(G)*∩D_(L)* where D_(P)*=c_(P) ⁻¹(1) is the design subspacethat is feasible with respect to pointwise constraint alone, and so on.The design problem is solved in two phases, depicted in FIG. 8:

Phase 1

Prune the design space from D to D_(P)*=c_(P) ⁻¹(1), i.e., solve thefollowing (simpler) problem:

-   -   Find D_(P)*⊆D, such that for all x∈Ω∈D_(P)*:

g _(i)(x)≤0, for 0<i≤n _(p).  (35)

According to embodiments described herein, the above problem can besolved by computing a maximal set Ω_(P)*:=maxD* in the partial orderingof designs via set containment.

Unfortunately, this is not possible for D_(G)* and D_(L)*. In mostcases, one can at best generate a finite sample of feasible designs thatare superior in some way.

Phase 2

Explore the pruned design space D_(P)* to find a sample D_(dom)*⊂D* of(locally) “Pareto-dominant” designs that satisfy the remainingconstraints:

-   -   Find D_(dom)*⊂D_(P)*, such that for all x∈Ω∈D_(opt)*;

g _(i)(Ω)≤0, for n _(P) <i≤n _(P) +n _(G),

g _(i)(x;Ω)≤0, for n _(P) +n _(G) <i≤n _(C),  (36)

Pareto-dominance of Ω∈D_(dom)* means that for some neighborhoodN(Ω)⊆D_(P)* and Ω∈N(Ω) in the pruned design space, within which no otherdesign is superior to Ω with respect to all objective functions ƒ₁, ƒ₂,. . . , ƒ_(n) _(O) :D→R. This can be posed as a minimization problemshown in (37).

$\begin{matrix}{{{Find}\mspace{14mu}\Omega} \in {D_{P}^{*}\mspace{14mu}{to}\mspace{14mu}\left\{ \begin{matrix}{{{minimize}\mspace{14mu}{f_{j}(\Omega)}},{{{for}\mspace{14mu} 0} < j \leq n_{O}},} \\{{subject}\mspace{14mu}{to}\mspace{14mu}{constraints}\mspace{14mu}{in}\mspace{14mu}{(36).}}\end{matrix} \right.}} & (37)\end{matrix}$

Here, N(Ω) is an open set in the induced Hausdorff topology of D_(P)*=

*(Ω*) (all solid subsets of the maximal element). The objectivefunctions define another partial ordering over the pruned design space,whose maximal elements are sought, i.e., D_(dom)*:=maxD.

The above problem may be solved by iterative optimization guided by TSF.The TSF is defined with respect to global objective functions ƒ_(j)(Ω)and global constraints g_(i)(Ω)≤0 and is penalized/filtered using localconstraints g_(i)(x; Ω)≤0. Since global optimization (N(Ω):=D_(P)*) isNP-hard, local optimality is used.

The example below shows how the design space can be pruned with respectto pointwise constraints (including set constraints) without prematureoptimization. The process is illustrated using examples fromkinematics-based constraints that are common in assembly, packaging, andmanufacturing.

The following results on the existence and uniqueness of maximalpointsets and their informational completeness as a representation forentire feasible design spaces are central to design space pruning.

Proposition 1 (Existence and Uniqueness) For every strictly local (i.e.,pointwise) constraint g_(i)*(x)≤0, its feasibility halfspace H_(i) has amaximal element Ω_(i)*=maxH_(i), defined implicitly by the following PMCtest:

$\begin{matrix}{{1_{\Omega_{i}^{*}}(x)}:=\left\{ \begin{matrix}1 & {{{{if}\mspace{14mu}{g_{i}^{*}(x)}} \leq 0},} \\0 & {{otherwise},}\end{matrix} \right.} & (38) \\{{i.e.},\mspace{14mu}{{\Omega_{i}^{*}:} = {\left\{ {x \in \Omega_{0}} \middle| {{g_{i}^{*}(x)} \leq 0} \right\}.}}} & (39)\end{matrix}$

The maximality is in terms of set containment, i.e., every satisfactorydesign is contained in the maximal element: Ω∈H_(i)⇒Ω⊆maxH_(i).

Proposition 2 (Completeness) For every strictly local (i.e., pointwise)constraint g_(i)*(

_(i)(x))≤0, its feasibility halfspace H_(i) contains every solidΩ⊆Ω_(i)*=maxH_(i), i.e., every solid subset of the maximal element isalso feasible: Ω⊆maxH_(i)⇒Ω∈H_(i).

In terms of predicates, the design subspace H_(i)=c_(i) ⁻¹(1) (whichsatisfies (16)) can now be represented by a single design Ω_(i)*=c_(i)⁻¹(1) (whose all points satisfy (18)). The maximal solid is thus acomplete representation of the feasibility halfspace as the collectionof all of its solid subsets denoted by H_(i)=

*(Ω_(i)*).

Here is an intuitive but simplified reasoning:

-   -   1. The set Ω_(i)* defined by (38) or (39) contains all points        that satisfy the constraints.    -   2. Every solid subset Ω⊆Ω_(i)* of the maximal set satisfies the        constraint, because all of its points satisfy the constraint        independently of the global shape.    -   3. Conversely, every feasible solid Ω∈H_(i) is the subset of        Ω_(i)*, because it only includes points that satisfy this        constraint independently of the global shape.

Note that the constraint's independence of the shape of Ω is crucial forthis to hold. For global or local constraints with dependency on Ωitself, attempting to write a PMC similar to (38) leads to a circulardefinition where the right-hand side depends on the set itself. Forexample, the (global or local) constraints on FEA and printabilityanalyses do not lead to maximal elements because their constraintsg_(i)(

_(i)(x; Ω))≤0 depend on particular design instances. It does make senseto define the maximal set of a feasible space (e.g., using theset-builder definition in (39)) in a way that it depends on a particularinstance of that space. On the other hand, the set constraints, such ascontainment under a prescribed motion, may give rise to maximalelements.

The above reasoning does not take topological regularization intoaccount—there is no reason for a maximal set obtained via (38) or (39)to be a solid, hence it may not itself be a valid design. However, theremay exist a valid maximal element obtained by regularizing (38) or (39).The correct definition is shown in (40).

Ω_(i) *:=ki{x∈Ω ₀ |g _(i)*(x)≤0}.  (40)

Here, k and i are the topological closure and interior operators,respectively.

Proposition 3 (Design Space Pruning) Given a number of pointwiseconstraints g_(i)*(x)≤0 for i=1, 2, . . . , n_(P), the feasible designspace D*, defined by intersecting all feasibility halfspaces H_(i)*, hasa maximal element Ω_(P)*=maxD_(P)* that satisfies the uniqueness andcompleteness properties, i.e., Ω⊆D_(P)*

Ω∈D_(P)*. It can be obtained by intersecting all maximal elementsΩ_(i)*=maxH_(i):

$\begin{matrix}{{\Omega_{P}^{*} = {\bigcap_{1 \leq i \leq n_{P}}\Omega_{i}^{*}}},{i.e.},{{1_{\Omega_{P}^{*}}(x)} = {\bigwedge\limits_{1 \leq i \leq n_{P}}{1_{\Omega_{i}^{*}}(x)}}},} & (41)\end{matrix}$

Here, the intersection/conjunction operators may need to be regularized.

To see why this is true, note that any query point's membership in an(unknown) feasible design can be tested against all m constraintsindependently of other points' membership. If g_(i)*(x)≤0 for i=1, 2, .. . , n_(P), the point can (but does not have to) be included in thedesign, for it to be feasible. The feasible design is hence a subset ofall points that satisfy all pointwise constraints.

The above result enables computing on design subspaces as first-classobjects. Computationally, the feasibility halfspaces are representeduniquely by their maximal elements. The pruning of halfspaces (abstractoperation) is implemented by intersecting maximal elements, i.e.,conjuncting point membership tests (i.e. testing whether each point inspace to see if it belongs to the design or not) defined by pointwiseconstraints (concrete algorithm) in an arbitrary order as seen in FIG.27.

According to embodiments described herein, design space pruning can beabstracted by intersecting the design space 1110 with feasibilityhalfspaces 1120, 1130 as shown in FIGS. 11A-11C. This is not acomputable operation in general. However, for pointwise constraints, itcan be computed by intersecting maximal elements in the physical space.

Next maximal elements for set constraints in are considered examplesfrom real-world engineering problems are provided.

Many design requirements, especially those that relate the interactionof moving shapes, may be expressed as set constraints of the formΓ(Ω)⊆E, where Γ:D→

(R³) is a set transformation. For example,

-   -   In packaging and assembly problems, the part's shape is often        designed so that it is restricted to remain within a specified        envelope while moving according to a prescribed motion.    -   When designing part to be machined using a given tool that moves        in the presence of obstacles (e.g., the part itself and        fixtures), the surface of the part has to be accessible without        collisions with the obstacles.

The common theme to many motion-based set transformation Γ:D→

(R³) is that they can refactored via (33) via a pointwise transformationγ:Ω₀→

(R³) that depends on the motion. The maximal pointset that satisfies thecontainment Γ(Ω)⊆E or non-interference test (Γ(Ω)∩E^(c))=∅ is definedimplicitly by its PMC given in terms of γ as shown in (43)-(45).

$\begin{matrix}{{1_{\Omega_{i}^{*}}(x)}:=\left\{ \begin{matrix}1 & {{{{if}\mspace{14mu}{\gamma(x)}} \subseteq E},} \\0 & {{otherwise},}\end{matrix} \right.} & (43) \\{{i.e.},\mspace{14mu}{{\Omega_{i}^{*}:} = \left\{ {x \in \Omega_{0}} \middle| {{\gamma(x)} \in E} \right\}}} & (44) \\{{= {{\gamma^{- 1}(E)}\bigcap\Omega_{0}}},} & (45)\end{matrix}$

Here, the set operators are regularized, as before. Note that γ⁻¹ maynot be a function, because γ is not necessarily invertible.

Here are a few classical examples from solid modeling:

-   -   For one-parametric sweep Γ(Ω):=sweep(M, Ω), one has γ(x)=Mx        where M=M(t)∈SE(3) is a continuous one-parametric set of motions        for t∈[t_(min), t_(max)]. The maximal shape that satisfies        Γ(x)⊆E is given by an unsweep unsweep(M, E).    -   For Minkowski sum Γ(Ω):=(Ω⊕B), one has γ(x)=(x+B) where B⊂R³ is        typically a solid. The maximal shape that satisfies Γ(x)⊆E is        given by a Minkowski difference with (E⊖(−B)).    -   For general dilation (which subsumes the above two) with general        rigid motions, the maximal shape is given by general erosion.    -   For non-rigid (but pointwise pre-determined) deformations, the        maximal shape is obtained by its PMC in terms of the pointwise        displacement function.

Procedure

Propositions 1 through 3 suggest a systematic procedure to prune thedesign space, i.e., reduce the design space D to D_(P)*=D∩(H₁∩H₂∩ . . .∩H_(n) _(P) ) for kinematic criteria expressed in terms of pointwise setconstraints:

-   -   Step 0. Initialize the feasible design space with the design        domain, i.e., in algorithmic terms, Ω_(P)*←Ω₀.    -   Step 1. Express the set constraint that outlines one of the        conditions for a given design Ω∈D to be feasible in the form        Γ(Ω)⊆E. Check if it can be restated as a pointwise constraint        γ(x)⊆E as in (32).    -   Step 2. Formalize the forward problem in terms of the PMC test        for the maximal element obtained from the pointwise constraint        as prescribed by (43).    -   Step 3. Invoke an IP-solver for the inverse problem, which        computes (an exact or approximate representation of) the maximal        element Ω_(i)* in (44).    -   Step 4. Prune the design space by intersecting the maximal        element Ω_(P)* obtained so far the new Ω_(i)*, i.e.,        Ω_(P)*←(Ω_(P)*∩Ω_(i)*) is a smaller maximal element representing        a pruned feasible subspace.    -   Repeat steps 1-4 for all pointwise set constraints.

Notice that the above procedure can be applied to different constraintsvia independent invocation of IP-solvers in an arbitrary orders.

The IP-solver in step 3 can be implemented in two fundamentallydifferent ways to obtain either an implicit representation (i.e., usingPMC test in (43)) or an explicit representation (i.e., using inversionin (45)) of the maximal element:

-   -   If there is access to an IP-solver that computes an explicit        representation (e.g., B-rep) of the inverse transformation        γ⁻¹(E), it can be directly intersected (using any CAD kernel)        with the design domain to obtain the maximal element as        prescribed by (45).    -   If there is access to a FP-solver that computes (explicitly or        implicitly) the forward transformation γ(x) for a given query        point x∈Ω₀, an approximate representation (e.g., point cloud or        voxelization) of the maximal element can be computed by:    -   (a) sampling the design domain with a sufficiently dense set of        query points;    -   (b) invoking the FP-solver to PMC-test them using (45); keep the        ones that pass the test and discard the ones that do not; and    -   (c) (optional) use adaptive local re-sampling around the points        that passed the test to obtain a better approximation.

Because of the independence of pointwise test, invocation of theFP-solver for different query points can be done with perfectparallelization.

The following process illustrates an example of pruning for containmentof moving parts in accordance with embodiments described herein.Consider the latch design problem introduced earlier, where the goal isto design a car hood latch that remains within an envelope E⊆Ω₀ whilemoving according to a motion M⊆SE(3).

Step 1

Every feasible latch design Ω∈D may satisfy the set constraint Γ(Ω)⊆Ewhere Γ(Ω):=sweep(M, Ω); i.e., the swept volume by the latch after beingtransformed by all configurations τ∈M (including any combination oftranslations and rotations, parametrized or otherwise) remains withinthe envelope. The sweep is a pointwise transformation, i.e., can becomputed as the union of all γ(x):=Mx=∪_(τ∈M) τx, which represents thetrajectory traced by the query point x∈Ω₀ along the prescribed motion.Hence, the containment constraint be tested in a pointwise fashion byγ(x)⊆E.

Step 2

Using the definitions in (43), a PMC test is constructed for the maximalshape in the design space that satisfies this pointwise constraint:

$\begin{matrix}{{1_{\Omega_{1}^{*}}(x)}:=\left\{ \begin{matrix}1 & {{{if}\mspace{14mu}{\forall{\tau \in {M:{{\tau\; x} \in E}}}}},} \\0 & {{otherwise}.}\end{matrix} \right.} & (46)\end{matrix}$

The forward problem involves following the trajectory, either exactly orapproximately (e.g., by sampling), and testing whether it remainsentirely within the envelope.

Step 3

The dual properties of the FP- and IP-solvers (i.e., Sweep and Unsweep)can be leveraged to construct an exact or approximate representation ofΩ₁*, as illustrated in FIGS. 12A-12D. An Unsweep solver may be used todirectly compute Ω₁*=unsweep(M⁻¹, E)∩Ω₀. However, if there is onlyaccess to efficient Sweep solver, an approximate representation of Ω₁*can be computed using the PMC test in (46) for a sufficiently densesample of query points and retaining the points whose forward trajectoryremains within the envelope as shown in (47) and (48).

[∀τ∈M:τx∈E]

x∈∩ _(τ∈M)τ⁻¹(E), i.e.,  (47)

sweep(M,Ω)⊆E

Ω⊆unsweep(M ⁻¹ ,E).  (48)

The invertibility of rigid transformations M

M⁻¹ may be key to efficient direct implementation of Unsweep.

FIGS. 12A-12D show unsweep(R, E) is the largest set that remains withinthe square envelope E 1210 while moving by a 21° clockwise rotation Raround a fixed pivot, i.e., without violating the containment constraintas depicted in 1220. The Unsweep can be computed in general byintersecting all moved instances of the envelope with an inverse motionas depicted in 1230. The result serves as a maximal shape 1240 forsubsequent material reduction (e.g., by TO).

Ω₁* can be sent as the initial design to design space exploration via TOor any other downstream material reducing IP-solver. It is guaranteedthat every valid design Ω⊆Ω₁* that is the output of the downstreamIP-solver, no matter how complicated, will continue to satisfy the setconstraint sweep(M, Ω)⊆E.

FIG. 13A shows a manufacturing setup with six clamps 1310, 1315, 1320,1325, 1330, 1340 used to locate and hold a designed part, and a 2-axisinstrument 1350 that can move in the plane. The envelope 1370 withinwhich the shape may be designed is shown. FIG. 13B illustrates theconvolution of the reflected head with the fixtures gives a 3D fieldwhose 2D cross-section via the plane of motion captures the collisionvolume at different translations. The position with zero collisionvolume 1380 are accessible. FIG. 13 C shows the 2D cross-section of the2.5D maximal manufacturable shape is obtained as the zero-set of theconvolution field. This pointset serves as the initial design for TO.

Consider another example in a very different setting. Suppose a part isbeing designed that may need to be fixtured in a crowded workholdingenvironment so that a machining instrument is able to access specificlocations without colliding with surrounding fixtures.

For simplicity, assume that the raw stock Ω₀ is a thick sheet of metal,fixtured on the machine bench. The manufacturing process is an EDMwire-cut or CO₂ laser-cut in which the tool assembly T=(H∪L) movesaccording to a planar (i.e., 2D) motion M⊆R² parallel to the workpiece.Here, it is assumed that the translation with a vector (x, y)∈M bringsthe wire or laser beam, abstracted by a vertical line of zero thicknessL, in contact with a line segment with x=(x, y, z)∈Ω₀ for a range ofz-values along the sheet thickness.

The head H cannot collide with the workpiece because they are located atdifferent elevations, i.e., (sweep(M, H)∩Ω)=∅ is a priori guaranteed,thus imposes no constraint. Nonetheless, the head H may collide with thefixtures F, which may extend above the workpiece, i.e., (sweep(M,H)∩F)=∅ imposes a constraint on the motion. This, in turn, imposes aconstraint for manufacturability, because the motion defines theboundary of the cut shape, restricted to a curve on the 2D plane:M=∂Ω∩R². In other words, designing the as-manufactured part's shapeamounts to designing the motion, because every translation of the wireor laser beam is in one-to-one correspondence with a point on the part'sboundary at which there is no collision between the head and thefixtures.

Step 1

The moving head H may not collide with the fixtures F when swept underthe motion M, i.e., (sweep(M, H)∩F)=∅, i.e., sweep(M, H)⊆F^(c). This iswritten in the standard form Γ(Ω)⊆E. Equivalently, for all translations(x, y)∈M, ((H+(x, y, 0))∩F)=∅, i.e., (H+(x, y, 0))⊆F^(c) which is inpointwise form γ(x)⊆E for x=(x, y, 0). Notice that γ(x): (H+x) does notdepend on Ω, as required.

Step 2

Using the definitions in (43), a PMC test is constructed for the maximalshape in the design space that satisfies this pointwise constraint:

1_(Ω) ₂ _(*)(x,y,z)=1 if(H+(x,y,0))⊆F ^(c),0 otherwise.  (49)

The success in defining a pointwise constraint (hence a PMC) depended onthe assumption of planar translation at a higher elevation than thepart, which guaranteed (H+(x, y, 0))⊆Ω^(c). Otherwise, the correctconstraint in (55) would have (F∪Ω)^(c)=(F^(c)∩Ω^(c)) instead of F^(c)on the right-hand side, making it stricter. But the maximal element Ω₂*and its corresponding feasibility halfspace H₂=

*(Ω₂*) cannot be defined in terms of a particular instance Ω∈H₂(circular definition).

Step 3

Once again, the dual properties of the FP- and IP-solvers can beleveraged to compute an explicit or implicit representation of the 2.5Dmaximal manufacturable solid.

The FP-solver can be any collision detection algorithm between arbitrarysolids, taking as input the displaced head (H+(x, y, 0)) above aparticular query point (x, y, z)∈Ω₀ and the stationary fixtures F.Therefore, one can sample the design domain (i.e., the raw stock) over a2D grid G and construct a bitmap image, representing the 2D section ofthe 2.5D maximal solid, by testing (55) for all (x, y)∈G. Every testrequires invoking the collision detection algorithm, and can be done inparallel.

Alternatively, one can construct an IP-solver to compute the collectionof all collision-free 2D translations (i.e., the configuration spaceobstacle). Since solids and regularized intersections exist, the setconstraint (H+(x, y, 0))∩F=∅ can be rewrittern in terms of measuresvol[(H+(x, y, 0))∩F]=∅ and convert it to an (in)equality constraint via(23). Hence, (49) becomes:

1_(Ω) ₂ _(*)(x,y,z)=¬sign∘vol[(H+(x,y,0))∩F]  (50)

=¬sign∫_(R) ₃ 1_(H)(x′−(x,y,0))1_(F)(x′)dv[x′],  (51)

Here, 1_(H+(x,y,0))(x′)=1_(H)(x′−(x, y, 0)) is the indicator function ofthe translated head. The integral on the right-hand side is aconvolution (1_(−H)*1_(F))(x), evaluated at x:=(x, y, 0), after areflection 1_(−H)(x′)=1_(H)(−x′). The integrand is nonzero only at x′∈R³where both indicator functions are nonzero, hence the integral does notvanish within the (measurable) regions of intersection. Substitutingthis relation into (57) yields:

$\begin{matrix}{{1_{\Omega_{2}^{*}}\left( {x,y,z} \right)} = \left\{ \begin{matrix}1 & {{{{if}\mspace{14mu}\left( {{\overset{\sim}{1}}_{H}*1_{F}} \right)\left( {x,y,0} \right)} = 0},} \\0 & {{otherwise}.}\end{matrix} \right.} & (52)\end{matrix}$

The convolution can be converted to pointwise multiplications in thefrequency domain using (forward and inverse) Fourier transforms:

(1_(−H)*1_(F))=

⁻¹{

{1_(−H)}·

{1_(F)}},  (53)

Here, (53) can be rapidly computed as fast Fourier transforms (FFT) andaccelerated on GPUs. FIG. 13A illustrates the results of thiscomputation. The convolution computes a 3D image in one shot using threeFFT computations on two 3D bitmaps (voxelized −H and F). However, onlyits 2D cross-section is may be needed at z=0, whose zero-set gives a 2Dbitmap image representation of the 2.5D maximal solid Ω₂*.

The remaining steps 4 and 5 are straightforward.

Note that everything in the above analysis would remain valid if themanufacturing instrument was allowed to rotate in the plane, except thatthe constraint on the right-hand side of (49) would have to be changedto hold for at least one planar rotation R∈SO(2) of the head H:

$\begin{matrix}{{1_{\Omega_{2}^{*}}(x)} = \left\{ \begin{matrix}1 & {{{{if}\mspace{14mu}{\exists{R \in {{SO}(2)}}}}:{\left( {{RH} + \left( {x,y,0} \right)} \right) \subseteq F^{c}}},} \\0 & {{otherwise}.}\end{matrix} \right.} & (54)\end{matrix}$

Accordingly, the convolution in (52) is adjusted:

$\begin{matrix}{{1_{\Omega_{2}^{*}}(x)} = \left\{ \begin{matrix}1 & {{{{{if}\mspace{14mu}{\exists{R \in {{SO}(2)}}}}:{\left( {1_{{- R}H}*1_{F}} \right)\left( {x,y,0} \right)}} = 0},} \\0 & {{otherwise}.}\end{matrix} \right.} & (55)\end{matrix}$

Here, the rotation can be parameterized as R=R(θ) for θ∈[0,2π) and1_(−RH)(x′)=1_(H)(−R⁻¹x′) where R⁻¹(θ)=R(−θ) is an inverse rotation. Tocompute the PMC, one has to sample the rotation angles θ∈[0,2π] and foreach trial rotation, resample the rotated head's 3D bitmap into the samegrid in which the fixtures' are rasterized to compute the discreteconvolution via FFT.

In more general manufacturing scenarios, a number of assumptions thatenabled pointwise formulation may be invalidated. For example, in a5-axis CNC machine, one deals with 6D rigid motions (R, t)∈SE(3)composed of 3D rotations R∈SO(3) and 3D translations t∈R³. The toolassembly T=(H∪C) is no longer guaranteed to avoid collisions with theworkpiece, leading to global constraints that depend on the part's shapeas well as the head, cutter, and fixtures. The configuration spaceobstacle M_(Ω):=obs(O_(Ω), T) where O_(Ω)=(Ω∪F) is stated as a groupconvolution ★ operation, which, in turn, can be computed as a Euclideanconvolution * as before as shown in (56).

1_(MΩ)(R,t)=sign∘(1_(O) _(Ω) ★1_(−T))(R,t)  (56)

=sign∘(1_(O) _(Ω) *1_(−RT))(t),  (57)

Attempting to write a PMC similar to (55) fails for several reasons; letus give it a try:

$\begin{matrix}{{1_{\Omega_{3}^{*}}(x)}\overset{?}{=}\left\{ \begin{matrix}1 & {{{{{if}\mspace{14mu}{\exists{R \in {{SO}(3)}}}}:{\left( {1_{O_{\Omega}}*1_{{- R}T}} \right)(t)}} = 0},} \\0 & {{otherwise}.}\end{matrix} \right.} & (58)\end{matrix}$

The first obvious problem is the dependency of the right-hand side on Ω,which makes for a circular definition. Moreover, the cutter's shapecannot be ignored (unlike the case with wire-/laser-cut). Hence, thereis no obvious way to assign a correspondence between the translationst∈R³ and the points x∈Ω₀ within the design domain, unless all possiblecontact configurations are considered and treat boundary pointsdifferently from interior points. Last but not least, passing acollision check at the contact configuration is not sufficient foraccessibility, because there may not exist a connected path from theinitial configuration of the tool assembly to the cutting pose ofinterest. For example, if a downstream TO creates cavities in the designin 3D, none of the will be accessible (unlike 2.5D).

Here, constraints that cannot be stated in pointwise form due to globaldependencies are examined. It is assumed that the design has been prunedfor all pointwise constraints to produce an initial designΩ_(P)*=(Ω₁*∩Ω₂*∩ . . . Ω_(n) _(P) *)⊆Ω₀ for design space explorationwith regards to the remaining (n−n_(P))=(n_(G)+n_(L)) global and/orlocal constraints.

In accordance with embodiments described herein, the phase 2 problem canbe solved by using design space exploration. The goal is not to proposenew optimization algorithms besides the many existing ones. Rather, ageneral strategy is proposed to deal with constraints that cannot bestated in a pointwise fashion, to guide gradient-descent optimization.

In order to move deterministically in the design space in directionsthat consistently reduce the violation of these constraints, theirsensitivities are quantified to hypothetical local changes in thedesign. Different gradient-like quantities can be defined for differentdesign representations. Here, it is demonstrated that the approachspecifically for defining, augmenting, and filtering topologicalsensitivity fields (TSF) with global and local constraints.

Fixed-point iteration (a.k.a Picard iteration) is an effective approachfor numerically solving multi-objective optimization problems, where theproblem is iteratively solved through series of outer- and inner-loops.As the value of each objective function is changed in the outer-loop,its value is kept fixed in the inner loop. The fixed objective functionsare treated as equality constraints for the single-objective inner-loopoptimization. Among the many popular approaches, a Pareto tracinglevelset TO approach (PareTO) is used, because it produces valid designs(i.e., solids) at all intermediate steps. It was shown in that Paretotracing can also be extended to density-based approaches such as solidisotropic material with penalization (SIMP).

For example, in classical TO, the goal is to obtain light-weight stiffstructures, leading to two competing objectives (mass and compliance)with a one-dimensional Pareto frontier, The problem can be formulated asfollows:

$\begin{matrix}{{Find}\left\{ \begin{matrix}{{{\Omega \in {D_{P}^{*}:{{minimize}\mspace{14mu}{\overset{¯}{V}}_{\Omega}\mspace{14mu}{and}\mspace{14mu} J_{\Omega}}}} = {\lbrack f\rbrack^{T}\left\lbrack u_{\Omega} \right\rbrack}},} \\{{{{subject}\mspace{14mu}{{{to}\mspace{14mu}\left\lbrack K_{\Omega} \right\rbrack}\left\lbrack u_{\Omega} \right\rbrack}} = \lbrack f\rbrack},}\end{matrix} \right.} & (59)\end{matrix}$

Here, the volume fraction V _(Ω):=vol[Ω]/vol[Ω_(P)*] is the ratio of the(unknown) vol[Ω] to the initial design's volume vol[Ω_(P)*], whereΩ_(P)*⊆Ω₀ is the maximal feasible pointset obtained from pruning.Classical TO in the absence of pruning is subsumed as a special casewhen Ω_(P)*=Ω₀. The second objective function J_(Ω)=[f]^(T) [u_(Ω)] isthe compliance (i.e., strain energy) obtained from FEA, in which [u_(Ω)]is the discretized displacement field and [f] is the external loadvector given as (Neumann) boundary conditions. The FEA also appears asan equality constraint [K_(Ω)][u_(Ω)]=[f] in which [K_(Ω)] is thestiffness matrix obtained from the design shape, material properties,and restraints given as (Dirichlet) boundary conditions.

The problem can be reformulated as a single-objective optimization for afixed volume fraction as seen in FIG. 28. Here, ILI stands forinner-loop iteration. Within each ILI, a single-objective optimizationis solved to minimize compliance J_(Ω) subject to a fixed volumefraction constraint V _(Ω)=V _(Ω) ^(targ) for a fixed 0<V _(Ω)^(targ)≤1. In PareTO, one starts off on the Pareto frontier at theright-most extreme with Ω:=Ω_(P)* and V _(Ω) ^(targ)=1, i.e., thebest-case scenario for compliance at the cost of the largest volume. Thealgorithm incrementally removes material to decrease V _(Ω) ^(targ) byintroducing holes in the design, without deviating too much from thePareto front. The ILI is a fixed-point iteration that applies localmodifications to the new design to bring it back to the Pareto front.

The inner-loop optimization can be expressed as local minimization ofthe augmented Lagrangian defined as:

_(Ω):=[f]^(T)[u _(Ω)]+λ₁( V _(Ω) −V _(Ω) ^(targ))+[λ₂]^(T)([K _(Ω)][u_(Ω)]−[f])  (61)

The Karush-Kuhn-Tucker (KKT) conditions [11] for this problem are givenby ∇

_(Ω)=0 in which the gradient is defined by partial differentiation withrespect to the independent variables; namely, the design variables usedto represent Ω and the Lagrange multipliers λ₁ and [λ₂]. The lattersimply encodes the constraints into ∇

_(Ω)=0:

$\begin{matrix}{{{\frac{\partial}{\partial\lambda_{1}}\mathcal{L}_{\Omega}} = {{\left( {{\overset{¯}{V}}_{\Omega} - {\overset{¯}{V}}_{\Omega}^{t\arg}} \right):} = 0}},} & (62)\end{matrix}$

$\begin{matrix}{{{\frac{\partial}{\partial\lambda_{2}}L_{\Omega}} = {{{{\left\lbrack K_{\Omega} \right\rbrack\left\lbrack u_{\Omega} \right\rbrack} - \lbrack f\rbrack}:} = \lbrack 0\rbrack}},} & \lbrack 63\rbrack\end{matrix}$

According to embodiments described herein, differentiation with respectto Ω∈D_(P)* depends on the particular parameterization used to representthe design by a finite set of decision variables for optimization. Thesevariables can be geometric/size variables (e.g., thickness in trussoptimization), density variables (e.g., volume fractions in SIMP), andso on. The goal is to present a representation-agnostic form in terms ofTSF.

A prime symbol (⋅)′ is used to represent the generic (linear)differentiation of a function with respect to Ω, (64) is obtained viathe chain rule.

$\begin{matrix}{{L_{\Omega\prime} = {{\lbrack f\rbrack^{T}\left\lbrack u_{\Omega\prime} \right\rbrack} + {\lambda_{1}{\overset{¯}{V}}_{\Omega\prime}} + {\left\lbrack \lambda_{2} \right\rbrack^{T}\left( {\left\lbrack K_{\Omega} \right\rbrack\left\lbrack u_{\Omega} \right\rbrack} \right)^{\prime}}}},{= {{\left( {\lbrack f\rbrack^{T} + {\left\lbrack \lambda_{2} \right\rbrack^{T}\left\lbrack K_{\Omega} \right\rbrack}} \right)\left\lbrack u_{\Omega\prime} \right\rbrack} + {\lambda_{1}{\overset{¯}{V}}_{\Omega\prime}} + {{\left\lbrack \lambda_{2} \right\rbrack^{T}\left\lbrack K_{\Omega\prime} \right\rbrack}\left\lbrack u_{\Omega} \right\rbrack}}},} & (64)\end{matrix}$

Computing [u_(Ω′)] is prohibitive, as it requires calling FEA as manytimes as the number of independent variables used to represent Ω. Thecommon solution is to choose [λ₂] such that [f]^(T)+[λ₂]^(T) [K_(Ω)]=[0](adjoint problem):

_(Ω′)=λ₁ V _(Ω′)+[λ₂]^(T)[K _(Ω′)][u _(Ω)], if[λ₂]:=−[K_(Ω)]⁻¹[f],  (65)

In general, if n_(O)>0 global objective functions ƒ_(j)(Ω) and anothern_(G)≥0 global (in)equality constraints g_(i)(Ω)≤0, (64) can begeneralized as:

_(Ω′):=Σ_(n) _(C) _(<i≤n) _(C) _(+n) _(O) λ _(i)ƒ_(i′)(Ω)+Σ_(n) _(P)_(<i≤n) _(P) _(+n) _(G) λ _(i) g _(i′)(Ω),   (66)

Here, the notation can be simplified by introducing F_(j)(Ω):=ƒ_(n) _(C)_(+j)(Ω)−ƒ_(n) _(C) _(+j) ^(targ) for 0<j≤n_(O) and F_(j)(Ω):=g_(n) _(P)_(−n) _(O) _(+j)(Ω) for n_(O)<j≤n_(O)+n_(G), hence:

_(Ω′)=Σ_(0<j≤n) _(O) _(+n) _(G) λ _(j) F _(j′)(Ω),  (67)

The ILI in FIG. 28 can be generalized to accommodate other globalobjective functions and global constraints. For example, suppose thepart is to be 3D printed at a given build direction. An additionalglobal (in)equality constraint is imposed in terms of an upper-bound V_(UB)≥0 on the total volume of support material that may be needed basedon an overhang angle criterion as seen in FIG. 29. Here, S_(≤⊆Ω) ^(c)represents the support structure. Its volume fraction V _(S) _(Ω)=vol[S_(Ω)]/vol[Ω_(P)*] can be computed as a function of the anglebetween surface normals and the build direction at every outer-loopiteration. The Lagrangian in (67) is further augmented by adding anotherterm λ₃(V _(S) _(Ω) −V _(UB)) and the generic sensitivity in (71) isupdated by incorporating V _(S) _(Ω) ′ as:

_(Ω′)=λ₁ V′ _(Ω)+[λ₂]^(T)[K′ _(Ω)][u _(Ω)]+λ₃ V _(S) _(Ω) ′.  (69)

FIG. 14 compares the solution designs based on a design domain 1410 to aTO problem with 1460 and without 1450 constraining the support materialvolume. Observe that optimization designs 1420, 1430, 1440 without thesupport constraint exits the feasibility halfspace with respect to thisconstraint for design volume fractions less than 70%. For lighterdesigns, the removed design material comes at the expense of additionalsupport material, hence costlier manufacturing. The fully constrainedoptimization designs 1425, 1435, 1445 with augmented sensitivity as in(69) dramatically increases the number of feasible and Pareto-optimaloptions, even at volume fractions lower than 70%. Here, FIG. 14 shows TOwith and without augmenting the sensitivity with constraints on thesupport material that may be needed for 3D printing along the verticalbuild direction. Many of the solutions without considering the supportconstraint will still satisfy that constraint due to the larger volumefraction occupied by the design itself. However, as the material isreduced from the design below 70%, the TO generates designs that requiremore support material.

Another example is TO subject to accessibility constraints formachining. Once again, a global (in)equality constraint can be imposedin (7) to in terms of an upper-bound V _(UB)≥0 on the total inaccessiblevolume as seen in FIG. 30. Here, R_(Ω)⊆Ω^(c) represents the maximalaccessible region outside the design for a combination of tools andapproach directions in 3-axis milling [50]. The volume fraction of theinaccessible regions is 1−V _(R) _(Ω) −V _(Ω) where V _(R) _(Ω)=vol[R_(Ω)]/vol[Ω_(P)*] and V _(Ω)=vol[Ω]/vol[Ω_(P)*], as before. TheLagrangian in (67) is further augmented by adding another term λ₃(1−V_(R) _(Ω) −V _(Ω)−V _(UB)), hence:

_(Ω′)=λ₁ V′ _(Ω)+[λ₂]^(T)[K′ _(Ω)][u _(Ω)]+λ₃( V′ _(R) _(Ω) −V′_(Ω))  (71)

=(λ₁−λ₃) V′ _(Ω)+[λ₂]^(T)[K′ _(Ω)][u _(Ω)]+λ₃ V′ _(R) _(Ω) ′.  (72)

One can alternatively formulate the optimization problem foraccessibility using the local constraint in as seen in FIG. 31. Here,the inaccessibility measure μ_(Ω)(x)=(1_(O) _(Ω) *{tilde over (1)}_(T))in (16), defined as the convolution in (17) is discretized to[μ_(Ω)]=[1_(O) _(Ω) *{tilde over (1)}_(T)] ands further simplified to[μ_(Ω)]=[1_(Ω)*{tilde over (1)}_(T)] assuming that the stationaryobstacle O_(Ω)=(Ω∪F) includes only the target design O_(Ω)=Ω, ignoringthe fixtures F:=∅. The tool assembly T=(H∪C) includes the holder H andcutter C, as before. Here, a conservative measure is used, aiming for noallowance for inaccessibility (i.e., μ₀:=0 in (16)) hence [1_(Ω)*{tildeover (1)}_(T)]=[0] over all discrete elements (e.g., voxels) whereverpossible in the design domain. The discrete convolution is computedusing two forward FFTs on [1_(Ω)] and [{tilde over (1)}_(T)], apointwise multiplication of their frequency domain grids, and an inverseFFT to obtain [1_(Ω)*{tilde over (1)}_(T)] in the physical domain (as avoxelized field).

Hereon, It is assumed that reducing the volume may be an objective/costfunction, hence the outer-loop is set up to incrementally decrease thevolume fraction budget V _(Ω) ^(targ)∈(0,1] starting from the initialvalue V _(Ω) ^(targ):=1 on the costlier extreme of the Pareto front. Theoptimization problem is formulated in general as shown in FIG. 32.

FIGS. 15A and 15B show that in the fixed-point iteration, the optimalityconditions are iteratively satisfied to ensure that at every step, thedesigns remain Pareto optimal. FIG. 15A shows an example Pareto tracingprocess in accordance with embodiments described herein. Starting fromthe initial design 1550, which may or may not be the maximal shapeobtained from design space pruning, the Pareto tracing approach removesmaterial from the least sensitive regions of the shape (as ranked byTSF) to obtain alternative Pareto-optimal designs 1540, 1530, 1525 alongthe Pareto front 1570. Every incremental step 1560 along the Paretofront in FIG. 15A involves a fixed-point iteration, illustrated in FIG.15B. The fixed-point iteration involves computing the Lagrangemultipliers and the TSF 1580, finding an updated iso-level set thatreduces the volume by a prescribed decrement 1582, and solving the FEAproblem on the updated shape 1584. The process is repeated until itconverges to a constant shape at the prescribed volume fraction.

TSF can be used to define

′_(Ω) and F′_(j)(Ω) in (67) in a representation-independent form. Let usfirst look at a few examples with manufacturability constraints inaddition to the physical constraints in FIG. 28.

The notion of a TSF is widely applied in the TO space as a means toguide the optimization process in moving from one candidate solution toanother in the search for local optima. Intuitively, the TSF is agradient-like operator for pointsets that quantifies the global effectof local changes of a given function (e.g., violation of a globalconstraint). The TSFs are coupled for various global and localconstraints in three distinct steps:

-   -   Defining TSFs for Global Constraints: For global constraints of        the general form g_(i)(Ω)≤0, one TSF per constraint is used to        measure how its violation changes after removing a hypothetical        small neighborhood (called an “inclusion”) at a given point.    -   Augmenting TSFs for Global Constraints: The individual TSFs are        linearly combined for all global constraints (including the        fixed objective functions).    -   Penalizing TSFs via Local Constraints: For local constraints of        the general form g_(i)(x; Ω)≤0, the TSF of the global        constraints is penalized by a linear combination of the        violation of local constraints.

For every function F_(j):D_(P)*→R that depends globally on the design(objective function or constraint), a field

: (Ω_(P)*×D_(P)*)→R is defined as its TSF via (75).

$\begin{matrix}{{{{\mathcal{T}_{j}\left( {x;\Omega} \right)}:} = {\lim\limits_{\epsilon\rightarrow 0^{+}}\frac{{F_{j}\left( {\Omega - {B_{\epsilon}(x)}} \right)} - {F_{j}(\Omega)}}{{vol}\left\lbrack {\Omega\bigcap{B_{\epsilon}(x)}} \right\rbrack}}},} & (75)\end{matrix}$

According to embodiments described herein, for 0<j≤n_(O)+n_(C).B_(∈)(x)⊂Ω₀ is a small 3D ball of radius ∈→0⁺ centered at a given querypoint x∈Ω. The numerator of the limit evaluates the (presumablyinfinitesimal) change in F_(j)(Ω) when the candidate design is modifiedas Ω

(Ω−B_(∈)(x)), i.e., by puncturing an infinitesimal cavity at the querypoint. The denominator vol[Ω∩B_(∈)(x)]=O(∈³) as ∈→0⁺ measures the volumeof the cavity. For internal points x∈iΩ (i.e., points that are notexactly on the topological boundary) one hasvol[Ω∩B_(∈)(x)]=vol[B_(∈)(x)] as ∈→0⁺.

The method of augmented Lagrangian can be extended to TSFs, and itseffectiveness was demonstrated by TO of multi-load structures underdeformation and stress constraints. The linear combination of thegeneric form in (67) is applied to compute an “augmented” TSF to couplethe global (in)equality constraints as shown in (76).

(x;Ω):=Σ_(0<j≤n) _(O) _(+n) _(G) λ_(j)

(x;Ω).  (76)

Note that the sum in (76) provides a representation-independentmathematical definition for the gradient in (67) with respect to the(unparameterized) pointset Ω∈Ω_(P)*. Rather than quantifying a directionof steepest descent for moving in a particular parameter space,

(x; Ω) identifies the set of points x∈Ω that are contributing the mostto the violation of constraints. A proper direction to move in the(unparameterized) design space D_(P)* is to remove the points withmaximal TSF.

The coefficients λ_(j)>0 has to be either computed by solving adjointproblems—as is shown for the case of strain energy in (65)—or selectedusing adaptive weighting schemes that are mainstream in multi-objectiveand multi-constraint TO.

The TSF operator maps global constraints to fields that vary dependingon x∈Ω. The local constraints are already defined as fields that vary ina similar fashion (i.e., are of the same “type” as the TSF). Hence, theTSF in (82) can be penalized with local constraints as:

(x;Ω):=

(x;Ω)+Σ_(n) _(P) _(+n) _(G) _(<i≤n) _(C) κ_(i) g _(i)(x;Ω).  (77)

The choice of coefficients κ_(i)>0 might need experimenting with the TOto adjust the relative importance of different constraints and improveconvergence properties.

The TSF orders the points in the design domain according to thepotential impact of removing their local neighborhoods on objectivefunction and constraints. An incremental improvement to the design isone that eliminates the points with the lowest TSF (e.g., the bottom5%). Here, ‘τ-modified’ (potentially infeasible) design Ω(τ)⊂Ω isdefined by a PMC in terms of the current design Ω:

$\begin{matrix}{{1_{\Omega{(\tau)}}(x)}:=\left\{ \begin{matrix}1 & {{{{if}\mspace{14mu}{\hat{\mathcal{T}}\left( {x;\Omega} \right)}} \geq \tau},} \\0 & {{otherwise},}\end{matrix} \right.} & (78) \\{{i.e.},\mspace{14mu}{{{\Omega(\tau)}:} = {\left\{ {x \in \Omega} \middle| {{\overset{\hat{}}{\mathcal{T}}\left( {x;\Omega} \right)} \geq \tau} \right\}.}}} & (79)\end{matrix}$

Here, the isolevel threshold τ>0 determines a step size for incrementalchange, e.g., τ:=0.05 means the least sensitive 5% are removed. It maybe important to select a small value so that only a small subset with

(x; Ω)≥τ is removed to obtain a shape that is not too different. The newdesign marginally violates the constraints and slightly deviates fromthe Pareto front 1580. However, it is close enough to the front that itcan be brought back by a fixed-point iteration 1582, 1584. The iterationmay not converge if the step size is too large. But if it does, itproduces another feasible and (locally) Pareto-dominant design that isslightly lighter.

Optimization Loops

Here is a general algorithm:

-   -   1. Pick a value δ>0 for the desired change in volume fraction        for the outer-loop iteration.    -   2. Compute        (x; Ω) from (76) and (77) and normalize it with its maximum        value over the current design.    -   3. Initialize Ω(τ)⊆Ω using (79) with a reasonably small initial        τ←τ₀ to start the fixed-point iteration:        -   (a) Cycle over the FP-solvers and update the performance            fields            _(i)(Ω)→            _(i)(Ω(τ)) (e.g., the constrained physical or kinematic            properties) for the τ-modified design obtained from (78).        -   (b) Re-evaluate the constraints using the updated            performance results; recompute the TSF using (76) and (77)            everywhere accordingly.        -   (c) Find τ>0 such that the τ-modified design in (86) with            the updated TSF has the desired reduction in volume            fraction, i.e., V _(Ω(τ))≈(V _(Ω)−δ).        -   (d) Repeat (a-c) until the τ-modified design does not            change. The result is feasible with respect to the            constraints and is Pareto-dominant.    -   4. Repeat (1-3) until the volume fraction reaches the smallest        feasible value to sustain the requirements.

Procedure

Here is a systematic procedure to explore the pruned design space, i.e.,trace a locally Pareto-optimal family of alternative design variantsD_(dom)*⊂D_(P)* by recurrent incremental thresholding of (augmented andpenalized) TSF, defined in terms of global and local constraints:

-   -   Step 0. Start at the extreme end of the Pareto front (maximal        volume) by initializing the design with the maximal pointset        obtained from pruning Ω←Ω_(P)*.    -   Step 1. Express the global objective functions and global and        local constraints for a given design Ω∈D_(P)* to formulate the        problem in the general form of (80).    -   Step 2. Define a subroutine to evaluate TSFs for each global        objective function and global constraint using (75), combine        them using (76), and penalize them with local constraints using        (77).    -   Step 3. Invoke the outer-loop optimization algorithm explained        above to incrementally reduce the material by thresholding the        TSF.    -   Step 4. Within the inner-loop (fixed point iteration) cycle over        FP-solvers to evaluate the objective functions and constraint        upon every incremental change in the outer-loop. Repeat until        the deviated solution converges back on the Pareto front.    -   Repeat steps 2-4 sequentially until the algorithm cannot find a        solution for after removing enough material, i.e., arrives at        other extreme end of the Pareto front (minimal volume).

Consider the car hood latch problem with the following kinematic andphysical constraints:

-   -   1. The latch may retain special features designated by the        designer, as illustrated in FIG. 16A.        -   One feature ensures that its mating pin (moving vertically            up and down) rotates the latch by 21° due to sliding contact            maintained through a spring (not shown here).        -   The other feature is for safety considerations; it ensures            that if the pin moves upwards in a sudden reverse motion—due            to a failure of the perimary latch—the secondary latch stops            claps it to prevent the car hood from opening.    -   2. As the latch rotates around its pivot from 0° to 21°, it may        remain completely within a safe region of space to avoid        interference with other car parts.    -   3. The latch is to be manufactured from stainless steel 304        using a metal AM process.    -   4. The latch should not weigh more than 0.30 pound.    -   5. The latch will experience loads at pre-determined        points/surfaces, including the contact forces with the pin        exerted by the spring. Under these loads, its maximum deflection        may not exceed 0.03 inches.

Such a diverse set of requirements is quite common, and should besimultaneously handled by the computational design framework. It isnoted from the first requirement that modeling design intent andsynthesizing functional features to satisfy them are difficult withoutknowing substantial information about the application. These featuresare given in a pre-processing step shown in FIG. 16A. Nevertheless, asubstantial remaining portion of the geometry is not defined byfunctional features and can be optimized. The remaining requirements aresystematically solved using embodiments described herein. FIG. 16A showsa preprocesses for a latch where functional surfaces are specified. FIG.16B illustrates that the Unsweep removes parts of the pre-processedinitial design that would exit the envelope for any clockwise rotationof θ∈[0°, 21°] around the pivot.

Step 0

As previously discussed, requirement 2 can be satisfied upfront (withoutpremature optimization) by pruning the design space via an IP-solver.The TO is started with the initial design Ω₌₁*:=unsweep(M, E) whereM={R(θ)∈SO(2)|0≤θ≤21°} is the collection of all rotations that the latchcan go experience, and E⊆R³ is the containment envelope (FIG. 13).

Since TO is a material-reducing procedure, the remaining requirements3-5 can be satisfied by TO without violating the containment constraint.

Steps 1, 2

In the absence of manufacturing constraints, the physics-basedconstraints for this problem are posed in common form of (59). Theupper-bound on the weight can be converted to an upper-bound on volumefraction V _(Ω)≤V _(Ω) ^(targ) where V _(Ω) ^(targ)=(0.30lb/ρ_(SS304))/vol[Ω₁*] using the known density of SS304.

The upper-bound on deflection _(Ω)(x)≤_(UB):=0.03″ need not be stated asa separate constraint, because it implies an upper-bound on compliance,hence a lower-bound on the volume fraction.

Steps 3, 4

At every outer-loop iteration, the maximal deflection increases due toremoved material. The algorithm checks if the deflection constraint isviolated and stops at the lightest possible solution.

Within the inner-loop fixed point iteration, the TSF is computed as in(65), based on which the τ-modified design Ω(τ) is extracted as theτ-superlevel set of the TSF. Subsequently, the FEA solver is invoked tosolve [K_(Ω(τ))][u_(Ω(τ))]=[f]. Based on the updated stiffness matrix[K_(Ω(τ))] and displacement field [u_(Ω(τ))] in response to the boundaryconditions, the Lagrange multipliers are updated via (65) as[λ₂]−[K_(Ω(τ)]) ⁻¹[f] and the TSF in (76) is recomputed. The iterationis repeated until the design remains unchanged.

FIG. 17 illustrates the Pareto front for solving the above problem,starting from the pruned design domain Ω_(P)* as prescribed above(strategy 1). The results of strategy 1 1750, 1755, 1760 are comparedagainst the the results 1730, 1735, 1740, 1745 where the algorithmstarts from the initial design domain Ω₀ (strategy 2), ignoring thecontainment constraint (requirement 2). The latter is an example ofpremature optimization, after which there is no guarantee that thedesign can be fixed to take the containment constraint into account. Thegraph also shows an incorrect attempt to fix the design as in FIG. 2,leading to an infeasible and suboptimal design (strategy 3) 1770. FIG.17 shows optimization fronts traced by different strategies. The firststrategy that applies TO to the pruned feasible design subspace of thecontainment constraint is the most computationally rational, as theentire family of solutions satisfy all requirements.

Let us next consider the problem of design for manufacturability withthe setup shown earlier. Once again, functional features are specifiedin pre-processing, as illustrated in FIG. 18. Here, it is assumed thatthese functional features are also accessible. In some cases, they couldbe introduced upfront in the raw stock and be modified by designerthrough trial and error. FIG. 19 illustrates the Pareto front ofaccessible designs optimized under specified loading boundary conditions1910 in accordance with embodiments described herein. The boundaryconditions 1910 are shown on the top-right corner of FIG. 19, includingboth forces and restrained surfaces. The underlying material isstainless steel with Young's modulus Y=200 GPa and Poisson's ratiov=0.33.

Step 0

The design space can be pruned with respect to the accessibility of a2-axis CNC instrument, where the 2D cross-section of the 2.5D maximalpointset was obtained as the 0-level set of a 3D convolution fieldbetween the head H and fixtures F (both in 3D), i.e.,Ω₂*(1_(−H)*1_(F))⁻¹(0).

Steps 1, 2

Here, the interest is finding a set of designs with maximal stiffnesswhile reducing the volume of the pruned design domain by another 60%.

All the optimized designs may have uniform cross-sections alongwire-/laser-cutting direction (i.e., are 2.5D). This constraint can beimposed either by applying a 2D TO to the cross-section of the initialdesign and extruding its results, or by using a 3D TO with a through-cutfiltering of the TSF. Here, the latter is used (PareTO in 3D) for thisexample.

The remaining steps 3 and 4 are similar to the previous example. Theonly FP-solver in the loop is a standard FEA, in absence of coupledmanufacturing constraints, noting that manufacturability is a prioriguaranteed in the pruning phase. In the next example, a problem isconsidered in which the manufacturing constraints cannot be pruned andhas to be coupled with the physical constraints within the inner-loopfixed-point iteration.

FIG. 20 shows the fixturing setup, raw stock, maximal manufacturabledomain (i.e., initial design for TO), and the optimized design at 40%volume fraction. FIG. 19 shows the Pareto front as it is traced from100% to 40% volume fraction. As with the previous example, the materialreducing nature of TO ensures that it does not violate themanufacturability constraint. FIG. 20 shows the optimized design at 40%volume fraction in the manufacturing setup of FIG. 13A.

Previously, two formulations for optimization subject to accessibilityconstraints were discussed, one with global formulation in FIG. 30(based on total inaccessible volume), and one with local formulation inFIG. 31 (based on inaccessibility measure as a convolution field).Experimenting with the former fails, as expected, because the TSF forthe global form is discontinuous. The convolution field, on the otherhand, is relatively well-behaved and can be used to penalize the TSF, asconfirmed by numerical experiments.

Step 0

In the absence of uncoupled pointwise constraints, start with theinitial design domain Ω₀∈D.

Step 1

The objective functions and constraints are given in (79). Theadditional constraint [1_(Ω)*{tilde over (1)}_(T)]=[0] requires thatevery point in the (voxelized) design [1_(Ω)] be accessible by the(voxelized) tool assembly [1_(T)], inverted as [{tilde over (1)}_(T)].Remember that the convolution's value at a given query point measuresthe volume of collision when the tool is displaced in such a way that arepresentative point on the tool—i.e., the origin of its localcoordinate system in which [1_(T)] is represented—is brought to thequery point in the design domain. The convolution field is defined overthe configuration space of relative motions (translations in this case).The proper selection of the local coordinate system is used to“register” the convolution field with the design domain and other fieldsdefined over it (e.g., TSF). The origin of the tool is picked at the tipof the cutter to simplify the formulation.

For the interior points, the constraint is violated, because the toolcannot reach the interior without colliding with the part. The violationis larger for points that are farther from the boundary, providing acontinuous penalty for TSF. Not every point in the exterior isaccessible either. Even if the tip of the cutter does not collide withthe part, the rest of the tool assembly might. The penalty is typicallysmaller for external points, as illustrated by FIGS. 21 and 25A-25C.when there are more than one tool or approach orientations, thealgorithm picks the minimum collision measure for penalization.

Step 2

The TSF for compliance is computed as usual, and is normalized by itsmaximum. An independent subroutine computes the convolution via FFTs, asdiscussed earlier, and is normalized by the volume of the tool(upper-bound to convolution). The TSF is penalized via convolution usingan adaptive weight (For instance start with λ₃: 0.01 and increase it toλ₃:=0.2 for lower volume fractions). Other design constraints such asminimum feature size or surface retainment can also be imposed.

Steps 3, 4

The outer-loop iteration is as before. The inner loop iteration nowcycles through one more FP-solver (the FFT-based convolution routine).

A simple 2D cantilever beam of FIG. 21 with simple boundary conditionsof a downward force F=1 N, Young's modulus of Y=1 GPa, and Poisson'sratio of v=0.3. Given a T-shaped tool with cutting part at the thin end,consider accessibility in two scenarios with the tool approaching fromone orientation (from left) and two orientations (from left and right).FIG. 21 illustrates how the convolution fields differs in the two cases.Subsequently, the compliance TSF is penalized to capture accessibilityunder tool orientations.

FIG. 21 shows penalizing TSF by the inaccessibility measure for T-shapedtool approaching from left 2110 or from both left and right 2120. Notethat the field is asymmetric for 2110. FIG. 22 illustrates optimizedshapes with and without the accessibility constraints with one 2220 andtwo 2230 tool orientations. Since the tool can only move in the plane,the TO cannot introduce interior holes without incurring a largepenalty. Moreover, it can only remove material from the boundary in sucha way that the remaining shape is machinable—e.g., no concave featuresof smaller feature size than the tool thickness in this case. Note thatthis is automatically enforced by penalizing convolution, withoutappealing explicitly to any notion of features or feature size. Here,optimized topologies at volume fractions 0.55 and 0.80 withoutaccessibility constraint 2210, with accessibility constraint for tool at0° 2220, and for tool at 0° and 180° 2230.

In the case of the tool approaching at 0°, material can be accessed andremoved only from the left side. However, with 0° and 1800 angles fortool orientation, the material can be removed from both sides. FIG. 23shows the Pareto fronts of the three scenarios 2310, 2320, 2330. Asexpected, optimization without accessibility constraints 2310 yields thebest performance in terms of compliance while imposing accessibilitywith one approach direction 2320 significantly increases the compliance.However, when the tool can approach from both directions 2330 comparableperformance to the unconstrained solutions can be achieved.

Consider the car hood latch example of FIGS. 16A and 16B. FIG. 24 showsthe optimized latches at 35% volume fraction with 2410 and without 2420the accessibility constraint. The same T-shaped tool is considered andoriented at both 0° and 180°. Imposing the accessibility constraintincreased the relative compliance from 1.09 to 1.26. FIG. 25A shows theoriginal TSF for compliance, FIG. 25B illustrates the inaccessibilitymeasure obtained from a convolution of the design and tool at 0° and180°, and FIG. 25C shows the penalized TSF to incorporate accessibilityand retain functional surfaces for the final design at a volume fractionof 35%.

It should be noted that in the above examples, only the collisionbetween the tool and the part is considered (i.e., no fixtures).Moreover, constraining the convolution field captures only the existenceof final collision-free configurations for the tool to machine the partat different points in the design domain. It does not guarantee acollision-free tool-path from the initial tool configuration to theremoval site.

For the definition of the TSF to be valid, the limit in (75) may existeverywhere in the design domain and for all intermediate designs. Inother words, puncturing the design with infinitesimal cavities may leadto infinitesimal changes in the violation of objective functions andconstraints. For this to hold, the functions may be sufficiently smooth.Here the function F_(j):D_(P)*→R may need to be differentiable in theHausdorff topology of D_(P)*=

*(Ω_(P)*), which is relative to the topology of D=

*(Ω₀). This is not always the case for general constraints. For example,if F_(j)(Ω) is itself evaluating a topological property, introducing apuncture (no matter how small) can produce a large change in F_(j)(Ω).For example, topological defects in AM due to resolution limits can becharacterized as integer-valued Euler characteristics. The Eulercharacteristic will increase by +1 after adding a cavity, hence

(x; Ω)˜+1/O(∈³)→+∞ as ∈→0⁺. In practice, this appears as a discontinuityin constraint evaluation which adversely affects the convergence of theoptimization loop (e.g., fixed-point iteration for PareTO).

As another example, recall the manufacturability constraint of (7) inwhich the total volume of inaccessible regions (for machining) wasupper-bounded as a global constraint. Puncturing a hole of volumevol[B_(∈)(x)] in the interior of the design adds exactly vol[B_(∈)(x)]to the inaccessible volume, hence

(x; Ω)=1 for all x∈iΩ. Most of the time, chipping off a visible (but toosmall) vol[Ω∩B_(∈)(x)] adds the same exact amount to the inaccessiblevolume, because the tool is of finite size and cannot remove that volumein practice. However, it is possible that a substantial region of designthat was initially inaccessible becomes accessible due to the smallchange at the boundary. Hence either

(x; Ω)=1 or

(x; Ω)→−∞ for all x∈∂Ω. This is not a well-behaved TSF for optimization.

The TSF is well-behaved for most global constraints that are defined asvolumetric integrals of continuous physical fields (e.g., strainenergy). This is the actual motivation behind using a volumetric measureof the inclusion in the denominator of (81) to normalize volumetricmeasures in the numerator. In principle, one can use a different measurefor global constraints that vary at a different rate than O(∈³). Smallchanges in the design often lead to small changes in physical response,when it is integrated over the entire domain, because integrationsmooths out the effects of local singularities. For instance, if aconstraint is imposed on the maximal stress as in (6), a realistic TSFcannot be defined due to stress concentration. The simplest model ofstress concentration yields an infinite maximal stressσ_(Ω)(x)˜O(∈^(−0.5)) near an infinitesimal radius of curvature ∈→0⁺,hence

(x; Ω)˜O(∈^(−3.5))→+∞ when defined for the maximal stress as in (6). Inpractice, one can alleviate this issue by using a volumetric integral(e.g., the p-norm) of stress rather than its maximum, in order to smoothout the effects of local singularities.

The other limitation with the usability of TSF is that the constraintfunction should not be locally “flat”. If the change in the violation ofa constraint is too small, i.e., decays faster than O(∈³), the limitwill vanish and the TSF will not help the optimization. For example,manufacturability constraints that have to do with surface propertiesare insensitive to volumetric changes in design.

In summary, the TSF formulation works well for global constraints thatchange smoothly with local volumetric material removal. Although this isnot true for all global constraints, the good news is that some of themcan be reformulated as a local constraints. Therefore, penalization ofthe local constraint can be used instead of defining a TSF for theglobal constraint. For example, although the accessibility constraintthat was discussed above does not yield a well-behaved TSF when treatedas a global constraint g₃(Ω)≤0 of (7), Previously, it was demonstratedthat it can be successfully incorporated to constraint PareTO usingpenalization of the local form g₄(x; Ω)≤0 of (16) in terms ofinaccessibility measure defined as the convolution in (17).

Mechanical design requires simultaneous reasoning aboutmultidisciplinary functional requirements and evaluating theirtrade-offs. These requirements are often expressed via heterogeneoustypes of constraints, including kinematics-based constraints forassembly and packaging, physics-based constraints for performance undermechanical or thermal loads, and both for manufacturability. Automateddesign optimization algorithms rarely consider all such requirements,and do not provide mechanisms to explore their trade space. For example,topology optimization can automatically generate designs with optimizedmaterial layouts for performance criteria such as strength andstiffness, but often ignores complex motion-based constraints imposed bycollision avoidance in assembly or accessibility in manufacturing.

The challenge in design space pruning and exploration is that IP-solversare usually equipped with the tools to satisfy only a subset of thecriteria in a multifunctional design problem. When these solvers arecomposed sequentially or in parallel, they can rarely provide guaranteesto retain the other criteria satisfied by the preceding or concurrentsolvers in the workflow. Moreover, most of the existing solvers generatea narrow subset of the design space-most commonly one or few design(s)that is/are deemed (locally or globally) “optimal” within the designsubspace that appears feasible to the solver. Such prematureoptimization dramatically limits the subsequent solvers' freedom toexplore (best-case scenario) and might even get deadlocked at infeasibledesigns (worst-case scenario).

To solve such challenges, a philosophy of treating design spaces (asopposed to individual designs) as first-class entities is followed—atleast to the extent that it is possible to do so by proper ordering ofsolvers in the workflow. This means that the entity being passed throughthe design pipeline—as input/output of consecutive synthesis solvers—isa design subspace described in its entirety by a representative object.This treatment allows postponing restrictive decisions and pushingpremature optimization downstream as much as possible. The designworkflows are organized by a careful analysis of the types ofconstraints and available solvers to address them and provide asystematic approach to compose FP- and IP-solvers depending on the typeof design constraint(s) they can satisfy.

A contribution of this work is a classification of constraints (namely,global, local, or strictly local) based on which the solvers areorganized systematically in the computational design workflow. Inparticular, the strictly local (i.e., pointwise) constraints can beevaluated without a knowledge of the global shape, hence lead to a pointmembership classification (PMC) for a maximal design that satisfiesthem. The maximal pointset represents the entire feasible designsubspace for a pointwise constraint, in the sense that containment inthe maximal pointset may be a necessary and sufficient condition forfeasibility. As such, the design space can be pruned upfront byintersecting maximal pointsets of pointwise constraints, withoutpremature optimization.

Most design criteria that depend on physics-based performance do notlead to a pointwise condition/PMC because the physical response of adesign at any given point is typically dependent on the overall shape,i.e., the membership of one point is couple with the membership of otherpoints. The dependency may be long-range, as in the case of staticequilibrium throughout a mechanical structure, or local, as in transientdynamic effects within a bounded neighborhood over a finite timeinterval. In either case, further design space pruning by means of PMC,to postpone decision making on the particular design layout, is not anoption. In such cases, the FP-solvers collaborate in generating feasibleand optimized designs by combining their sensitivity fields and methodssuch as fixed-point iteration for tracing the trade space of multipleobjectives.

A revelation of the classification is that the two different types ofproblems (namely, design space ‘pruning’ and ‘exploration’) demonstratea different form of duality between forward and inverse problem(IP/FP)-solvers for generative design:

-   -   For pointwise constraints, an IP-solver can be constructed from        an FP-solver by generating a large sample of points in the        design domain, applying the FP-solver in a pointwise fashion,        evaluating the constraint, and retaining/discarding the ones the        do/do not satisfy the constraint. In other words, the FP-solver        provides a PMC test for the IP-solver. The process can be        perfectly parallelized.    -   For other constraints, an IP-solver can be constructed from an        FP-solver by generating a number of candidate designs,        evaluating the constraints, obtaining a sensitivity field to        order the different points in the design according to their        expected impact on the (dis)satisfaction of the constraint,        removing the least sensitive points, and try again with the        modified design. In other words, the FP-solver provides ae        evaluator for candidate designs to put in a feedback loop for        the generate-and-test IP-solver. The process is a sequential        loop that is repeated until convergence.

A limitation of the approach is that it does not provide any guaranteesfor satisfying constraints that are neither pointwise nordifferentiable. For some local constraints (e.g., accessibility measuresfor machining), it has been shown that penalizing the sensitivity fieldsof other global constraints with the local constraint can be effective.

However, it is unclear how to systematically make such decisions withevery new problem and constraint, unlike the case with pointwiseconstraints (pruned upfront) or differentiable non-pointwise constraints(filtered via local sensitivity analysis).

The above-described methods can be implemented on a computer usingwell-known computer processors, memory units, storage devices, computersoftware, and other components. A high-level block diagram of such acomputer is illustrated in FIG. 26. Computer 2600 contains a processor2610, which controls the overall operation of the computer 2600 byexecuting computer program instructions which define such operation. Thecomputer program instructions may be stored in a storage device 2620(e.g., magnetic disk) and loaded into memory 2630 when execution of thecomputer program instructions is desired. Thus, the steps of the methodsdescribed herein may be defined by the computer program instructionsstored in the memory 2630 and controlled by the processor 2610 executingthe computer program instructions. The computer 2600 may include one ormore network interfaces 2650 for communicating with other devices via anetwork. The computer 2600 also includes a user interface 2660 thatenables user interaction with the computer 2600. The user interface 2660may include I/O devices 2662 (e.g., keyboard, mouse, speakers, buttons,etc.) to allow the user to interact with the computer. Such input/outputdevices 2662 may be used in conjunction with a set of computer programsas an annotation tool to annotate images in accordance with embodimentsdescribed herein. The user interface also includes a display 2664 fordisplaying images and spatial realism maps to the user. According tovarious embodiments, FIG. 26 is a high-level representation of possiblecomponents of a computer for illustrative purposes and the computer maycontain other components.

Unless otherwise indicated, all numbers expressing feature sizes,amounts, and physical properties used in the specification and claimsare to be understood as being modified in all instances by the term“about.” Accordingly, unless indicated to the contrary, the numericalparameters set forth in the foregoing specification and attached claimsare approximations that can vary depending upon the desired propertiessought to be obtained by those skilled in the art utilizing theteachings disclosed herein. The use of numerical ranges by endpointsincludes all numbers within that range (e.g. 1 to 5 includes 1, 1.5, 2,2.75, 3, 3.80, 4, and 5) and any range within that range.

The various embodiments described above may be implemented usingcircuitry and/or software modules that interact to provide particularresults. One of skill in the computing arts can readily implement suchdescribed functionality, either at a modular level or as a whole, usingknowledge generally known in the art. For example, the flowchartsillustrated herein may be used to create computer-readableinstructions/code for execution by a processor. Such instructions may bestored on a computer-readable medium and transferred to the processorfor execution as is known in the art. The structures and proceduresshown above are only a representative example of embodiments that can beused to facilitate embodiments described above.

The foregoing description of the example embodiments have been presentedfor the purposes of illustration and description. It is not intended tobe exhaustive or to limit the inventive concepts to the precise formdisclosed. Many modifications and variations are possible in light ofthe above teachings. Any or all features of the disclosed embodimentscan be applied individually or in any combination, not meant to belimiting but purely illustrative. It is intended that the scope belimited by the claims appended herein and not with the detaileddescription.

What is claimed is:
 1. A method comprising: receiving manufacturingcriteria for a product part; sorting the manufacturing criteria intodifferent classes of one or both objective functions and constraintsbased on when they can be satisfied or optimized; determining constraintviolations; and producing a design workflow to generate one or moredesigns of a part to comply with one or more of satisfying constraintsand optimizing objective functions, the design workflow invoking acombination of solvers including forward problem solvers and inverseproblem solvers.
 2. The method of claim 1, wherein the constraintscomprise one or more of set constraints, equality constraints, andinequality constraints.
 3. The method of claim 1, wherein scopes of theconstraints comprise one or more of global, local, and strictly local.4. The method of claim 1, further comprising computing performancefields of the one or more designs.
 5. The method of claim 4, furthercomprising: evaluating the objective functions based on the performancefields; and ordering the one or more designs based on the evaluatedobjective functions.
 6. The method of claim 4, further comprising:computing constraint violations based on the performance fields; anddetermining if the one or more designs are feasible based on thecomputed constraint violations.
 7. The method of claim 1, wherein thesolvers are organized in the workflow such that the solvers that producea broadest family of designs are invoked first.
 8. A method comprising:receiving performance criteria for a product part; sorting theperformance criteria into different classes of one or both objectivefunctions and constraints based on when they can be satisfied oroptimized; determining constraint violations; and producing a designworkflow to generate one or more designs of a part to comply with one ormore of satisfying constraints and optimizing objective functions, thedesign workflow invoking a combination of solvers including forwardproblem solvers and inverse problem solvers.
 9. The method of claim 8,wherein the constraints comprise one or more of set constraints,equality constraints, and inequality constraints.
 10. The method ofclaim 8, wherein scopes of the constraints comprise one or more ofglobal, local, and strictly local.
 11. The method of claim 8, furthercomprising computing one or more performance fields of the one or moredesigns.
 12. The method of claim 11, further comprising: evaluating theobjective functions based on the performance fields; and ordering theone or more designs based on the evaluated objective functions.
 13. Themethod of claim 11, further comprising: computing constraint violationsbased on the performance fields; and determining if the one or moredesigns are feasible based on the computed constraint violations. 14.The method of claim 8, wherein the solvers are organized in the workflowsuch that the solvers that produce a broadest family of designs areinvoked first.